大型数组或列表的 4 桶直方图的微观优化

And*_*eas 1 c# optimization simd histogram micro-optimization

我有一个特别的问题。我将尝试尽可能准确地描述这一点。

我正在做一个非常重要的“微优化”。一次运行数天的循环。所以如果我能减少这个循环时间,它需要一半的时间。10 天将减少到只有 5 天等。

我现在拥有的循环是函数:“testbenchmark1”。

我有 4 个索引需要在这样的循环中增加。但是当从列表中访问索引时,实际上需要一些额外的时间,正如我所注意到的。这就是我想知道是否有其他解决方案。

indexes[n]++; //increase correct index

“testbenchmark1”的完整代码需要 122 毫秒:

void testbenchmark00()
{
    Random random = new Random();
    List<int> indexers = new List<int>();
    for (int i = 0; i < 9256408; i++)
    {
        indexers.Add(random.Next(0, 4));
    }
    int[] valueLIST = indexers.ToArray();


    Stopwatch stopWatch = new Stopwatch();
    stopWatch.Start();

    int[] indexes = { 0, 0, 0, 0 };
    foreach (int n in valueLIST) //Takes 122 ms
    {
        indexes[n]++; //increase correct index
    }

    stopWatch.Stop();
    MessageBox.Show("stopWatch: " + stopWatch.ElapsedMilliseconds.ToString() + " milliseconds");
}
Run Code Online (Sandbox Code Playgroud)

现在下面的“testbenchmark2”代码只是实验性的,我知道它不正确,但我想知道是否有任何类似的方法来使用这种数字:“1_00_00_00_00”,如果有可能看到:“00_00_00_00”为四个不同的整数。例如,如果我要进行求和:1_00_00_00_00 + 1_00_01_00_00 = 1_00_01_00_00然后最后可以提取每个数字,四个中的每一个都像这样:00, 01, 00, 00

但我不知道即使使用二进制数是否有可能。是的任何类型的解决方案。只需添加这样的数字。就像测试一样,循环只用了 59 毫秒,是 122 毫秒时间的一半。所以我很有趣,看看是否有任何想法?

double num3 = 1_00_00_00_00;
double num4 = 1_00_01_00_00;
for (int i = 0; i < valueLIST.Count; i++) //Takes 59 ms
{
    num3 += num4;
}
Run Code Online (Sandbox Code Playgroud)

“testbenchmark2”的完整代码需要 59 毫秒:

void testbenchmark2()
{
    List<String> valueLIST = new List<String>(); 
    for (int i = 0; i < 9256408; i++) //56
    {
        valueLIST.Add(i.ToString());
    }

    //https://www.geeksforgeeks.org/binary-literals-and-digit-separators-in-c-sharp/
    double num3 = 1_00_00_00_00;
    double num4 = 1_00_01_00_00;

    Stopwatch stopWatch = new Stopwatch();
    stopWatch.Start();
    for (int i = 0; i < valueLIST.Count; i++) //Takes 59 ms
    {
        num3 += num4;
    }
    stopWatch.Stop();
    MessageBox.Show("stopWatch: " + stopWatch.ElapsedMilliseconds.ToString() + " milliseconds\n\n" + num3);
}
Run Code Online (Sandbox Code Playgroud)

编辑
下面是我正在尝试做的更干净的代码!
但是下面的代码可能是正确的或解决方案,但它显示了我相信我尝试做的事情。

        void newtest()
        {
            double num1 = 1_00_00_00_00;
            double num2 = 1_00_01_00_00;
            double num3 = 1_00_01_01_00;

            List<double> testnumbers = new List<double>();
            testnumbers.Add(num1);
            testnumbers.Add(num2);
            testnumbers.Add(num3);

            double SUM = 0;
            for (int i = 0; i < testnumbers.Count; i++)
            {
                SUM += testnumbers[i];
            }

            //The result is
            //300020100

            //Would it possible to extract the "four buckets" that I am interesting in somehow?
            //00_02_01_00
        }
Run Code Online (Sandbox Code Playgroud)

Pet*_*des 6

在使用 AVX2 的现代 x86-64(如 Skylake 或 Zen 2)上,每 2.5 个时钟周期(每个内核)大约 8 个元素(1 个 AVX2 向量)应该是可能的。或每 2 个时钟展开。或者在您的 Piledriver CPU 上,使用 AVX1 可能每 3 个时钟有 1 个 16 字节的索引向量_mm_cmpeq_epi32

一般策略适用于 2 到 8 个存储桶。对于字节、16 位或 32 位元素。(因此,字节元素为您提供每 2 个时钟周期直方图的 32 个元素,最佳情况是在溢出之前收集字节计数器的一些外循环开销。)

更新:或者将一个 int 映射到1UL << (array[i]*8)使用 SIMD / SWAR 添加来增加计数器的 4 个字节之一,我们可以接近 1 个时钟在 SKL 上的每个 8 个 int 向量,或在 Zen2 上每 2 个时钟。(这更特定于 4 个或更少的存储桶和 int 输入,并且不会缩小到 SSE2。它需要变量移位或至少 AVX1 变量洗牌。)将字节元素与第一种策略一起使用可能仍然更好就每个周期的元素而言。

正如@JonasH 指出的那样,您可以让不同的内核处理输入数组的不同部分。在典型的台式机上,单核可以接近饱和内存带宽,但众核至强具有较低的每核内存带宽和更高的聚合,并且需要更多内核来饱和 L3 或 DRAM 带宽。 为什么 Skylake 在单线程内存吞吐量方面比 Broadwell-E 好得多?


一次运行数天的循环。

在迭代非常慢的单个输入列表上,它仍然不会溢出 int 计数器?或者重复调用不同的大列表(比如你的 ~900k 测试数组)?

我相信我想避免增加列表或数组的索引,因为它似乎消耗了很多时间?

这可能是因为您在禁用优化的情况下进行基准测试。不要那样做,一点意义都没有;通过禁用优化,不同的代码会减慢不同的速度。更明确的步骤和 tmp 变量通常会使调试模式代码变慢,因为需要使用调试器查看更多内容。但是当您使用正常优化进行编译时,它们可以优化为正常的指针增量循环。

遍历数组可以有效地编译成 asm。

缓慢的部分是通过内存的依赖链,用于增加数组的变量索引。例如,在 Skylake CPU 上,add具有相同地址的内存目标以每 6 个时钟周期大约一个增量重复瓶颈,因为下一个add必须等待加载前一个存储的值。(来自存储缓冲区的存储转发意味着它不必先等待它提交到缓存,但它仍然比添加到寄存器慢得多。)另请参阅 Agner Fog 的优化指南:https ://agner.org /优化/

由于计数仅分布在 4 个存储桶中,您将在很多情况下等待重新加载另一条最近指令存储的数据,因此如果计数良好,您甚至无法实现每个时钟周期接近 1 个元素分布在 L1d 缓存中仍然很热的更多计数器上。

这个问题的一个很好的解决方案是用多个计数器数组 展开循环在 SIMD 中矢量化直方图的方法?. 就像int[] indexes = { 0, 0, 0, 0 };您可以将其设为每个包含四个计数器的 2D 数组一样。您必须手动展开源中的循环以遍历输入数组,并处理展开部分之后剩下的最后 0..3 个元素。

对于中小型计数数组来说,这是一种很好的技术,但如果复制计数器开始导致缓存未命中,则该技术会变得很糟糕。


使用窄整数来节省缓存占用空间/内存带宽。

您可以/应该做的另一件事是为 0..3 个值的数组使用尽可能窄的类型:每个数字都可以放在一个字节中,因此使用 8 位整数将为您节省 4 倍的缓存占用空间/内存带宽.

x86 可以有效地向/从完整寄存器加载/存储字节。使用 SSE4.1,您还可以使用 SIMDpmovzxbd来提高在循环中byte_array[i]使用a 时自动矢量化的效率int_array[i]

(当我说 x86 时,我的意思是包括 x86-64,而不是 ARM 或 PowerPC。当然,您实际上并不想编译 32 位代码,Microsoft 称之为“x86”)


桶的数量非常少,比如 4

这看起来像是 SIMD 比较的工作。使用 x86 SSE2,int每 16 字节数据向量的元素数等于您的直方图箱数。

您已经有了 SIMD 类型的想法,试图将数字视为四个单独的字节元素。请参阅https://en.wikipedia.org/wiki/SIMD#Software

00_01_10_11它只是数字中人类可读分隔符的源级语法,并且double是一种浮点类型,其内部表示与整数不同。而且你绝对不想使用字符串;SIMD 可以让你同时操作一个整数数组的 4 个元素。

我能看到的最好方法是分别为 4 个值中的每个值计算匹配项,而不是将元素映射到计数器。 我们希望并行处理多个元素,但是当一个元素向量中存在重复值时,将它们映射到计数器可能会发生冲突。您需要将该计数器增加两次。

等价的标量是:

int counts[4] = {0,0,0,0};
for () {
    counts[0] += (arr[i] == 0);
    counts[1] += (arr[i] == 1);
    counts[2] += (arr[i] == 2);  // count matches
  //counts[3] += (arr[i] == 3);  // we assume any that aren't 0..2 are this
}
counts[3] = size - counts[0] - counts[1] - counts[2];
// calculate count 3 from other counts
Run Code Online (Sandbox Code Playgroud)

这(在 C++ 中)GCC-O3实际上会按照我在下面手动执行的方式自动矢量化https : //godbolt.org/z/UJfzuH。Clang 甚至在自动矢量化时展开它,所以它应该比我的手动矢量化输入版本更好int。尽管如此,仍然不如vpermilps这种情况的替代策略好。

(如果您想要具有有效窄总和的字节元素,您仍然需要手动矢量化,仅在外循环中加宽。)


对于字节元素,请参阅如何使用 SIMD 计算字符出现次数。元素尺寸对于计数器来说太窄了;它会在 256 次计数后溢出。因此,您必须在内循环中加宽,或者使用嵌套循环在加宽之前进行一些累积。

我不懂 C#,所以我可以用 x86 程序集或带有内在函数的 C++ 编写代码。也许 C++ 内在函数对您更有用。C# 有某种向量扩展,应该可以移植它。

这是用于 x86-64 的 C++,使用 AVX2 SIMD 内在函数。有关一些信息,请参阅https://stackoverflow.com/tags/sse/info

// Manually vectorized for AVX2, for int element size
// Going nearly 4x as fast should be possible for byte element size

#include <immintrin.h>

void count_elements_avx2(const std::vector<int> &input,  unsigned output_counts[4])
{
    __m256i  counts[4] = { _mm256_setzero_si256() };  // 4 vectors of zeroed counters
                  // each vector holds counts for one bucket, to be hsummed at the end

    size_t size = input.size();
    for(size_t i = 0 ; i<size ; i+=8) {  // 8x 32-bit elements per vector
        __m256i v = _mm256_loadu_si256((const __m256i*)&input[i]);  // unaligned load of 8 ints
        for (int val = 0 ; val < 3; val++) {
           // C++ compilers will unroll this with 3 vector constants and no memory access
            __m256i match = _mm256_cmpeq_epi32(v, _mm256_set1_epi32(val));  // 0 or all-ones aka -1
            counts[val] = _mm256_sub_epi32(counts[val], match);   // x -= -1 or 0 conditional increment
        }
    }


    // transpose and sum 4 vectors of 8 elements down to 1 vector of 4 elements
    __m128i summed_counts = hsum_xpose(counts);   // helper function defined in Godbolt link
    _mm_storeu_si128((__m128i*)output_counts, summed_counts);

    output_counts[3] = size - output_counts[0]
                       - output_counts[1] - output_counts[2];

    // TODO: handle the last size%8 input elements; scalar would be easy
}
Run Code Online (Sandbox Code Playgroud)

这与 clang 编译得很好(在Godbolt 编译器资源管理器上)。据推测,您可以编写编译为类似机器代码的 C#。如果不是,请考虑从 C++ 编译器调用本机代码(如果无法从编译器获得真正最佳的代码,则使用 asm 手写)。如果您的实际用例运行的迭代次数与您的基准测试一样多,那么如果不必复制输入数组,则可以分摊额外的开销。

 # from an earlier version of the C++, doing all 4 compares in the inner loop
 # clang -O3 -march=skylake
.LBB0_2:                                     # do {
    vmovdqu ymm7, ymmword ptr [rcx + 4*rdx]    # v = load arr[i + 0..7]
    vpcmpeqd        ymm8, ymm7, ymm3           # compare v == 0
    vpsubd  ymm4, ymm4, ymm8                   # total0 -= cmp_result
    vpcmpeqd        ymm8, ymm7, ymm5
    vpsubd  ymm2, ymm2, ymm8
    vpcmpeqd        ymm7, ymm7, ymm6           # compare v == 2
    vpsubd  ymm1, ymm1, ymm7                   # total2 -= cmp_result
    add     rdx, 8                             # i += 8
    cmp     rdx, rax
    jb      .LBB0_2                          # }while(i < size)
Run Code Online (Sandbox Code Playgroud)

估计的最佳 Skylake 性能:每个向量约 2.5 个周期(8 int 或 32 int8_t)

或 2 与展开。

没有 AVX2,只使用 SSE2,你会有一些额外的movdqa指令,每个向量只能做 4 个元素。不过,这仍然是内存中的标量直方图的胜利。即使是 1 个元素/时钟也不错,并且应该可以使用可以在任何 x86-64 CPU 上运行的 SSE2。

假设没有缓存未命中,硬件预取到 L1d 保持在循环之前。这可能只会发生在 L2 缓存中已经很热的数据上。 我还假设内存对齐没有停顿;理想情况下,您的数据按 32 字节对齐。 如果通常不是,如果数组足够大,可能值得处理第一个未对齐的部分,然后使用对齐的负载。

对于字节的元素,最内环路将类似于(具有vpcmpeqbvpsubb但hsumming到64位计数器,以避免溢出之前仅运行最多255(未256)迭代。因此可以通过每个矢量将是相同的,但与4倍每个向量尽可能多的元素。

有关性能分析的详细信息,请参阅https://agner.org/optimize/https://uops.info/。例如vpcmpeqd在 uops.info

Haswell/Skylake 的内部循环只有 9 个融合域 uops,因此最好的前端瓶颈是每 2.25 个周期约 1 次迭代(管道宽度为 4 个 uops)。 小循环效应在某种程度上受到了影响:在执行 uop 计数不是处理器宽度倍数的循环时,性能是否会降低?- Skylake 的循环缓冲区因错误的微代码更新而被禁用,但即使在此之前,9 uop 循环最终平均每 2.25 个周期发出比 1 个迭代略差的结果,比方说 2.5 个周期。

Skylakevpsubd在端口 0、1 或 5 上运行,并vpcmpeqd在端口 0 或 1 上运行。因此端口 0、1、5 上的后端瓶颈是 3 个端口的 6 个向量 ALU uop,或每 2 个周期 1 次迭代。 所以前端瓶颈占主导地位。 (即使没有展开,Ice Lake 更宽的前端也可能让它在后端出现瓶颈;除非您使用 AVX512,否则后端吞吐量相同......)

如果 clang 从数组的末尾开始索引并将索引向上计数到零(因为它无论如何都选择使用索引寻址模式)它可以节省一个 uop,总共 8 个 uops = 每 2 个周期一次迭代在前面-end,匹配后端瓶颈。(无论哪种方式,标量add和宏融合cmp/jcc,或add/jcc循环分支都可以在端口 6 上运行,并且负载不会与 ALU 端口竞争。)即使在缓存未命中时,依赖于负载的 ALU uops 的 Uop 重放也不应该成为问题,如果 ALU uops 是瓶颈,通常会有大量旧的 uops 等待执行单元准备就绪,而不是等待加载数据。

展开 2 将具有相同的好处:分摊 2 uop 的循环开销。所以 2 个输入向量有 16 个 uops。 这是 SKL 和 IceLake 上的管道宽度以及 Zen 上的单 uop 管道宽度的好倍数。展开更多可以让前端保持在执行之前,但即使是任何后端延迟也会让前端在调度程序中建立一个 uop 缓冲。这将让它足够早地执行加载。

Zen2 具有更宽的前端(6 uop 或 5 条指令宽,IIUC)。这些指令都不是多 uop,因为 Zen2 将向量 ALU 扩展到 256 位,所以这是 5 个单 uop 指令。 vpcmpeq*在 FP 0,1 或 3 上运行,与 相同vpsubd,因此后端瓶颈与 Skylake 相同:每 2 个周期 1 个向量。但是更广泛的前端消除了这个瓶颈,即使没有展开,关键路径仍然是后端。

Zen1 每个 256 位向量操作需要 2 uop(或更多用于车道交叉,但这些是简单的 2 uop)。所以大概 12/3 = 每个包含 8 或 32 个元素的向量有 4 个周期,假设它可以有效地通过前端获得这些 uops。

我假设通过计数向量的 1 周期延迟依赖链由后端很好地安排并且不会导致许多浪费的周期。可能没什么大不了的,特别是如果您在现实生活中遇到任何内存瓶颈。(在 Piledriver 上,SIMD 整数操作有 2 个周期的延迟,但是可以运行它们的 2 个向量 ALU 端口的 6 个 ALU uops 是每 3 个周期 1 个向量(128 位),因此即使没有展开,也有足够的工作来隐藏该延迟。)

我没有分析这个的横向总和部分。它在循环之外,因此每次调用只需运行一次。您确实标记了这个微优化,但我们可能不需要担心那部分。


其他桶数

该策略的基本情况是 2 个桶:计数一件事的匹配项,count_other = 大小 - 计数。

我们知道每个元素都是这 4 种可能性之一,因此我们可以假设任何x不是 0、1 或 2 的元素都是3,而无需检查。这意味着我们不必指望比赛3在所有的,并且可以得到计数从桶size - sum(counts[0..2])

(在进行此优化之前,请参阅上述性能分析的编辑历史记录。我在进行此优化并更新 Godbolt 链接后更改了数字,希望我没有遗漏任何内容。)


Skylake-Xeon 上的 AVX512

对于 64 字节向量,无法vpcmpeqd创建全零 (0) 或全一 (-1) 元素的向量。相反,您将与掩码寄存器进行比较,并使用它来对set1(1). 喜欢c = _mm512_mask_add_epi32(c, _mm512_set1_epi32(1))

不幸的是,对比较结果位掩码进行标量 popcount 效率不高。


随机代码审查:在你的第一个基准测试中:

int[] valueLIST = indexers.ToArray();

这似乎毫无意义;根据 MS 的文档(https://docs.microsoft.com/en-us/dotnet/standard/collections/),列表可以有效地索引。我认为它相当于 C++ std::vector<T>。您可以只迭代它而无需复制到数组。


Alt 策略 - 将 0..3 映射到 int 的一个字节中的一组位

如果您不能将元素缩小到字节以节省内存带宽,那很好。

但是说到这里,在使用 3x pcmpeqb / psubb 进行计数之前,使用 2x _mm256_packs_epi32(vpackssdw) 和_mm256_packs_epi16( vpacksswb) 缩小到 8 位整数可能是值得的。每 4 个输入向量需要 3 个 uops,以将字节元素压缩为 1。

但是,如果您的输入以 int 元素开头,则最好不要打包然后比较 3 种方式。

您有 4 个存储桶,一个int有 4 个字节。 如果我们可以将每个int元素转换1为适当字节底部的a ,那么_mm256_add_epi8在扩展到 64 位计数器之前,我们最多可以添加255 次内循环迭代。(使用_mm256_sad_epu8对零技巧的标准对无溢出的无符号字节进行hsum。)

有两种方法可以做到这一点。第一种:使用 shuffle 作为查找表。AVX2vpermd工作 ( _mm256_permutexvar_epi32) 使用数据作为索引向量和一个常量_mm256_set_epi32(0,0,0,0, 1UL<<24, 1UL<<16, 1UL<<8, 1UL<<0)作为被打乱的数据。或者输入双关语向量以将 AVX1vpermilps用作 LUT,而 LUT 向量的上半部分也包含这些字节。

vpermilps更好:AMD Zen 1 上的 uops 更少,而且所有地方的延迟都更低,因为它是 in-lane。(可能会导致某些 CPU 的旁路延迟,从而减少延迟优势,但仍不比 差vpermd)。

出于某种原因vpermilps,矢量控制在 Zen2 上有 2 个周期的吞吐量,即使它仍然是单个 uop。或 Zen1 上的 4 个周期(对于 2 uop YMM 版本)。在英特尔上是 1 个周期。 vpermdAMD 的情况更糟:更多的 uops 和同样糟糕的吞吐量。

vpermilps xmm根据 Agner Fog 的测试,Piledriver 上的(16 字节向量)具有 1/clock 的吞吐量,并在“ivec”域中运行。(因此,当在“预期的”浮点操作数上使用时,它实际上具有额外的旁路延迟延迟,但在整数上则没有)。

   // Or for Piledriver, __m128 version of this

    __m256 bytepatterns = _mm256_casts256_ps(_mm256_set_epi32(
         1<<24, 1<<16, 1<<8, 1<<0,
         1<<24, 1<<16, 1<<8, 1<<0) );
    __m256i v = _mm256_loadu_si256((const __m256i*)&input[i]);
    v = _mm256_castps_si256(_mm256_permutevar_ps(bytepatterns, v));  // vpermilps 32-bit variable shuffle
    counts = _mm256_add_epi8(counts, v);

     // after some inner iterations, separate out the 
     // set1_epi32(0x000000ff) counts, 0x0000ff00 counts, etc.
Run Code Online (Sandbox Code Playgroud)

这将在每个int元素内产生交错计数器。如果您没有在 256 次计数之前累积它们,它们就会溢出。有关使用单个计数器的简单版本,请参阅如何使用 SIMD 计算字符出现次数

在这里,我们可能会展开并使用 2 个不同的 LUT 向量,因此当我们想要将所有计数0组合在一起时,我们可以2 个向量混合在一起并屏蔽掉其他向量。


除了改组,我们还可以使用 AVX2 变量转换来实现。

sums += 1UL << (array[i]*8); 其中*8是一个字节中的位数,也是用移位完成的。我把它写成一个标量 C++ 表达式,因为现在你有机会看到你的整数字节的想法是如何真正起作用的。只要我们不让单个字节溢出,SIMD 字节是在字节之间添加块进位还是使用 32 位双字元素都没有关系。

我们将使用 AVX2 这样做:

__m256i v = loadu...();
v = _mm256_slli_epi32(v, 3);  // v *= 8
v = _mm256_sllv_epi32(_mm256_set1_epi32(1), v);
counts = _mm256_add_epi8(counts, v);
Run Code Online (Sandbox Code Playgroud)

这是 2 个移位指令加上vpaddb. 在 Skylake 上,可变计数班次vpsllvd很便宜:单 uop 并在多个端口上运行。但在 Haswell 和 Zen 上,速度较慢。(与vpermilpsAMD相同的吞吐量)

对于 shuffle 版本,2 个端口的 2 uop 仍然无法击败 1 个端口的 1 uop。(除非您交替使用这两种策略来将工作分配到SKL上的所有 ALU 端口。)

因此,无论哪种方式,最内层循环都可以每个时钟运行 1 个向量,或者通过仔细交错移位与洗牌方法可能会稍微好一些。

但它需要在 128 或 255 次内部循环迭代中分摊少量开销。

最后的清理可能会将 2 个向量混合在一起以获得一个仅包含 2 个桶的计数的向量,然后vpshufb( _mm256_shuffle_epi8) 将同一桶的字节计数器分组到相同的 qwords 中。然后vpsadbw( _mm256_sad_epu8) 对零可以水平求和每个 qword for 中的那些字节元素_mm256_add_epi64。所以外循环工作应该是 2 vpblendw, 2x vpshufb, 2x vpsadbw, 2xvpaddq然后回到内循环的另外 255 次迭代。可能还会检查您是否在数组末尾的 255 次迭代内,以设置内部迭代的循环边界。