展开循环以实现Ivy Bridge和Haswell处理器的最大吞吐量

17

我正在使用AVX一次计算八个点积。在我的当前代码中,我做的类似于以下操作(未展开之前):

Ivy-Bridge/Sandy-Bridge

__m256 areg0 = _mm256_set1_ps(a[m]);
for(int i=0; i<n; i++) {        
    __m256 breg0 = _mm256_load_ps(&b[8*i]);
    tmp0 = _mm256_add_ps(_mm256_mul_ps(arge0,breg0), tmp0); 
}

哈斯韦尔

__m256 areg0 = _mm256_set1_ps(a[m]);
for(int i=0; i<n; i++) {      
    __m256 breg0 = _mm256_load_ps(&b[8*i]);
    tmp0 = _mm256_fmadd_ps(arge0, breg0, tmp0);
}

每种情况下需要展开循环多少次才能保证最大吞吐量?

对于使用FMA3的Haswell,我认为答案在这里 FLOPS per cycle for sandy-bridge and haswell SSE2/AVX/AVX2。 我需要展开循环10次。

对于Ivy Bridge,我认为是8。 这是我的逻辑。 AVX加法的延迟为3,乘法的延迟为5。 Ivy Bridge可以同时使用不同的端口进行一个AVX乘法和一个AVX加法。 使用符号m表示乘法,a表示加法,x表示未操作以及数字来表示局部和(例如,m5表示与第5个局部和相乘),我可以写成:

port0:  m1  m2  m3  m4  m5  m6  m7  m8  m1  m2  m3  m4  m5  ... 
port1:   x   x   x   x   x  a1  a2  a3  a4  a5  a6  a7  a8  ...

通过在九个时钟周期后使用8个部分和(四个来自加载和五个来自乘法),我可以每个时钟周期提交一个AVX加载、一个AVX加法和一个AVX乘法。

我想这意味着在Ivy Bridge和Haswell的32位模式中无法实现此任务的最大吞吐量,因为32位模式只有八个AVX寄存器?

编辑:关于悬赏问题。 我主要的问题仍然存在。 我想获得上面提到的Ivy Bridge或Haswell函数的最大吞吐量,n可以是大于或等于64的任何值。 我认为这只能通过展开(Ivy Bridge展开8次,Haswell展开10次)来完成。 如果您认为可以用另一种方法来解决,请让我们看看。 在某种程度上,这是如何实现每个周期的理论最大4 FLOP的变体?。 但是,我不仅要求乘法和加法,还要求每个时钟周期具有一个256位加载(或两个128位加载)、一个AVX乘法和一个AVX加法,或者使用Haswell每个时钟周期具有两个256位加载和两个FMA3指令。

我还想知道需要多少寄存器。 对于Ivy Bridge,我认为是10个。 一个用于广播,一个用于加载(由于寄存器重命名只有一个),八个用于八个部分和。 因此,我认为这不能在32位模式下完成(实际上,当我在32位模式下运行时,性能会显著下降)。

我应该指出,编译器可能会给出误导性的结果高度优化的矩阵乘法代码的MSVC和GCC之间性能差异

我目前在Ivy Bridge上使用的函数如下。 这基本上将一个64x64矩阵a的一行与所有64x64矩阵b相乘(我在a的每一行上运行此函数64次,以获取矩阵c中的完整矩阵乘积)。

#include <immintrin.h>
extern "C" void row_m64x64(const float *a, const float *b, float *c) {      
    const int vec_size = 8;
    const int n = 64;
    __m256 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
    tmp0 = _mm256_loadu_ps(&c[0*vec_size]);
    tmp1 = _mm256_loadu_ps(&c[1*vec_size]);
    tmp2 = _mm256_loadu_ps(&c[2*vec_size]);
    tmp3 = _mm256_loadu_ps(&c[3*vec_size]);
    tmp4 = _mm256_loadu_ps(&c[4*vec_size]);
    tmp5 = _mm256_loadu_ps(&c[5*vec_size]);
    tmp6 = _mm256_loadu_ps(&c[6*vec_size]);
    tmp7 = _mm256_loadu_ps(&c[7*vec_size]);

    for(int i=0; i<n; i++) {
        __m256 areg0 = _mm256_set1_ps(a[i]);

        __m256 breg0 = _mm256_loadu_ps(&b[vec_size*(8*i + 0)]);
        tmp0 = _mm256_add_ps(_mm256_mul_ps(areg0,breg0), tmp0);    
        __m256 breg1 = _mm256_loadu_ps(&b[vec_size*(8*i + 1)]);
        tmp1 = _mm256_add_ps(_mm256_mul_ps(areg0,breg1), tmp1);
        __m256 breg2 = _mm256_loadu_ps(&b[vec_size*(8*i + 2)]);
        tmp2 = _mm256_add_ps(_mm256_mul_ps(areg0,breg2), tmp2);    
        __m256 breg3 = _mm256_loadu_ps(&b[vec_size*(8*i + 3)]);
        tmp3 = _mm256_add_ps(_mm256_mul_ps(areg0,breg3), tmp3);   
        __m256 breg4 = _mm256_loadu_ps(&b[vec_size*(8*i + 4)]);
        tmp4 = _mm256_add_ps(_mm256_mul_ps(areg0,breg4), tmp4);    
        __m256 breg5 = _mm256_loadu_ps(&b[vec_size*(8*i + 5)]);
        tmp5 = _mm256_add_ps(_mm256_mul_ps(areg0,breg5), tmp5);    
        __m256 breg6 = _mm256_loadu_ps(&b[vec_size*(8*i + 6)]);
        tmp6 = _mm256_add_ps(_mm256_mul_ps(areg0,breg6), tmp6);    
        __m256 breg7 = _mm256_loadu_ps(&b[vec_size*(8*i + 7)]);
        tmp7 = _mm256_add_ps(_mm256_mul_ps(areg0,breg7), tmp7);    
    }
    _mm256_storeu_ps(&c[0*vec_size], tmp0);
    _mm256_storeu_ps(&c[1*vec_size], tmp1);
    _mm256_storeu_ps(&c[2*vec_size], tmp2);
    _mm256_storeu_ps(&c[3*vec_size], tmp3);
    _mm256_storeu_ps(&c[4*vec_size], tmp4);
    _mm256_storeu_ps(&c[5*vec_size], tmp5);
    _mm256_storeu_ps(&c[6*vec_size], tmp6);
    _mm256_storeu_ps(&c[7*vec_size], tmp7);
}

1
请注意,在Core i7及其后续版本中,展开小循环可能对性能产生负面影响。在Nehalem上,循环流探测器只能缓存28个µop,我不确定在Ivy Bridge/Haswell中这个尺寸是否增加了。 - Paul R
2
你使用展开技术的目的是什么?是为了避免易于预测的分支吗?虽然OOO能够跨多个迭代执行,但你仍应该只限于tmp0依赖。 - Leeor
4
他所说的并不是uop缓存,而是循环流探测器,它更小(且在带宽方面更有效)。 - Leeor
2
我也在质疑展开像FMAs这样延迟很长的依赖链的必要性。静态准备好索引的好处微不足道,你的性能受向量单元带宽/延迟的控制 - 计算索引(和预测一些分支)可以并行完成而不会产生影响。只有通过分析才能确定谁是正确的。 - Leeor
2
那只适用于Sandy/Ivy Bridge。Haswell可以每个周期执行两个256位的加载。我是在指你最初的循环,其中(我认为)你为每个FMA都有一个广播和一个加载。 - Mysticial
显示剩余23条评论
2个回答

16

对于Sandy/Ivy Bridge,您需要按3展开循环:

  • 只有FP Add依赖于循环的前一次迭代
  • FP Add可以每个周期发出指令
  • FP Add需要三个周期才能完成
  • 因此,按3展开循环/1 = 3将完全隐藏延迟
  • FP Mul和FP Load不依赖于前一次迭代,您可以依靠OoO核心以接近最佳顺序发出它们。这些指令可能会影响展开因子,但仅在它们降低FP Add的吞吐量时才会影响(在这里不是这种情况,FP Load + FP Add + FP Mul可以每个周期发出)。

对于Haswell,您需要按10展开循环:

  • 只有FMA依赖于循环的前一次迭代
  • FMA可以每个周期双重发出(即独立指令平均需要0.5个周期)
  • FMA的延迟为5
  • 因此,按5 / 0.5 = 10展开循环将完全隐藏FMA的延迟
  • 两个FP Load微操作不依赖于前一次迭代,并且可以与2x FMA一起协同发出,因此它们不影响展开因子。

我明白逻辑。 我不清楚OoO逻辑如何解决这个问题,但我知道在代码中该怎么做。 我将不得不阅读英特尔手册以了解更多信息。 由于我已经安排好了内存和64x64矩阵的方式,所以展开8次可能会产生最佳结果。 32位速度较慢,因为我需要超过八个寄存器。 但这仍然意味着,在Haswell上,无法在32位模式下获得最大吞吐量。 - Z boson
我刚意识到我的问题中实际上没有展开循环。我正在执行八个8宽度的点积(64个点积)。我使用八个AVX累加器读取/写入64侧矩阵的一行。虽然如此,这并不真正改变问题。问题只是变成了需要多少个累加器才能实现最大吞吐量,答案是相同的。 - Z boson
为什么Haswell上的mulss只需要3个周期,而不同于Agner指令表中的数据?(使用多个累加器展开FP循环)表明使用比最低要求更多的累加器有所帮助。显然,调度并不总是完美的,和/或者使用更多的累加器可以更好地吸收内存延迟变化。 - Peter Cordes

7

我只是在这里回答自己的问题以添加信息。

我已经对Ivy Bridge代码进行了分析。当我最初在MSVC2012中测试时,展开超过两个并没有帮助太多。然而,根据我在高度优化的矩阵乘法代码的MSVC和GCC之间性能差异中的观察,我怀疑MSVC没有最优地实现内部函数。因此,我使用g++ -c -mavx -O3 -mabi=ms在GCC中编译内核,将对象转换为COFF64并将其放入MSVC中。我现在得到的结论是,将展开次数增加到三次可以获得最佳结果,这证实了Marat Dunkhan的答案。

下面是在Xeon E5 1620 @3.6GHz MSVC2012中以秒为单位的时间。

unroll    time default            time with GCC kernel
     1    3.7                     3.2
     2    1.8 (2.0x faster)       1.6 (2.0x faster)
     3    1.6 (2.3x faster)       1.2 (2.7x faster)
     4    1.6 (2.3x faster)       1.2 (2.7x faster)

以下是在运用GCC和fma指令时,i5-4250U在Linux下的性能表现时间(g++ -mavx -mfma -fopenmp -O3 main.cpp kernel_fma.cpp -o sum_fma):

unroll    time
     1    20.3
     2    10.2 (2.0x faster)
     3     6.7 (3.0x faster) 
     4     5.2 (4.0x faster)
     8     2.9 (7.0x faster)
    10     2.6 (7.8x faster)

下面的代码适用于Sandy-Bridge/Ivy Bridge。对于Haswell,请使用例如tmp0 = _mm256_fmadd_ps(a8,b8_1,tmp0)
kernel.cpp
#include <immintrin.h>

extern "C" void foo_unroll1(const int n, const float *b, float *c) {      
    __m256 tmp0 = _mm256_set1_ps(0.0f);
    __m256 a8 = _mm256_set1_ps(1.0f);
    for(int i=0; i<n; i+=8) {
        __m256 b8 = _mm256_loadu_ps(&b[i + 0]);
        tmp0 = _mm256_add_ps(_mm256_mul_ps(a8,b8), tmp0);
    }
    _mm256_storeu_ps(c, tmp0);
}

extern "C" void foo_unroll2(const int n, const float *b, float *c) {
    __m256 tmp0 = _mm256_set1_ps(0.0f);
    __m256 tmp1 = _mm256_set1_ps(0.0f);
    __m256 a8 = _mm256_set1_ps(1.0f);
    for(int i=0; i<n; i+=16) {
        __m256 b8_1 = _mm256_loadu_ps(&b[i + 0]);
        tmp0 = _mm256_add_ps(_mm256_mul_ps(a8,b8_1), tmp0);
        __m256 b8_2 = _mm256_loadu_ps(&b[i + 8]);
        tmp1 = _mm256_add_ps(_mm256_mul_ps(a8,b8_2), tmp1);
    }
    tmp0 = _mm256_add_ps(tmp0,tmp1);
    _mm256_storeu_ps(c, tmp0);
}

extern "C" void foo_unroll3(const int n, const float *b, float *c) { 
    __m256 tmp0 = _mm256_set1_ps(0.0f);
    __m256 tmp1 = _mm256_set1_ps(0.0f);
    __m256 tmp2 = _mm256_set1_ps(0.0f);
    __m256 a8 = _mm256_set1_ps(1.0f);
    for(int i=0; i<n; i+=24) {
        __m256 b8_1 = _mm256_loadu_ps(&b[i + 0]);
        tmp0 = _mm256_add_ps(_mm256_mul_ps(a8,b8_1), tmp0);
        __m256 b8_2 = _mm256_loadu_ps(&b[i + 8]);
        tmp1 = _mm256_add_ps(_mm256_mul_ps(a8,b8_2), tmp1);
        __m256 b8_3 = _mm256_loadu_ps(&b[i + 16]);
        tmp2 = _mm256_add_ps(_mm256_mul_ps(a8,b8_3), tmp2);
    }
    tmp0 = _mm256_add_ps(tmp0,_mm256_add_ps(tmp1,tmp2));
    _mm256_storeu_ps(c, tmp0);
}

extern "C" void foo_unroll4(const int n, const float *b, float *c) {      
    __m256 tmp0 = _mm256_set1_ps(0.0f);
    __m256 tmp1 = _mm256_set1_ps(0.0f);
    __m256 tmp2 = _mm256_set1_ps(0.0f);
    __m256 tmp3 = _mm256_set1_ps(0.0f);
    __m256 a8 = _mm256_set1_ps(1.0f);
    for(int i=0; i<n; i+=32) {
        __m256 b8_1 = _mm256_loadu_ps(&b[i + 0]);
        tmp0 = _mm256_add_ps(_mm256_mul_ps(a8,b8_1), tmp0);
        __m256 b8_2 = _mm256_loadu_ps(&b[i + 8]);
        tmp1 = _mm256_add_ps(_mm256_mul_ps(a8,b8_2), tmp1);
        __m256 b8_3 = _mm256_loadu_ps(&b[i + 16]);
        tmp2 = _mm256_add_ps(_mm256_mul_ps(a8,b8_3), tmp2);
        __m256 b8_4 = _mm256_loadu_ps(&b[i + 24]);
        tmp3 = _mm256_add_ps(_mm256_mul_ps(a8,b8_4), tmp3);
    }
    tmp0 = _mm256_add_ps(_mm256_add_ps(tmp0,tmp1),_mm256_add_ps(tmp2,tmp3));
    _mm256_storeu_ps(c, tmp0);
}

main.cpp

#include <stdio.h>
#include <omp.h>
#include <immintrin.h>

extern "C" void foo_unroll1(const int n, const float *b, float *c);
extern "C" void foo_unroll2(const int n, const float *b, float *c);
extern "C" void foo_unroll3(const int n, const float *b, float *c);
extern "C" void foo_unroll4(const int n, const float *b, float *c);

int main() {
    const int n = 3*1<<10;
    const int r = 10000000;
    double dtime;
    float *b = (float*)_mm_malloc(sizeof(float)*n, 64);
    float *c = (float*)_mm_malloc(8, 64);
    for(int i=0; i<n; i++) b[i] = 1.0f;

    __m256 out;
    dtime = omp_get_wtime();    
    for(int i=0; i<r; i++) foo_unroll1(n, b, c);
    dtime = omp_get_wtime() - dtime;
    printf("%f, ", dtime); for(int i=0; i<8; i++) printf("%f ", c[i]); printf("\n");

    dtime = omp_get_wtime();    
    for(int i=0; i<r; i++) foo_unroll2(n, b, c);
    dtime = omp_get_wtime() - dtime;
    printf("%f, ", dtime); for(int i=0; i<8; i++) printf("%f ", c[i]); printf("\n");

    dtime = omp_get_wtime();    
    for(int i=0; i<r; i++) foo_unroll3(n, b, c);
    dtime = omp_get_wtime() - dtime;
    printf("%f, ", dtime); for(int i=0; i<8; i++) printf("%f ", c[i]); printf("\n");

    dtime = omp_get_wtime();    
    for(int i=0; i<r; i++) foo_unroll4(n, b, c);
    dtime = omp_get_wtime() - dtime;
    printf("%f, ", dtime); for(int i=0; i<8; i++) printf("%f ", c[i]); printf("\n");
}

我不确定这是否有帮助,但您可以在参数b和c上添加__restrict。 - user678269
乘以1f有什么作用? - user678269
为什么在只使用一次的情况下,b8_1...b8_4 这些表达式要使用变量? - user678269
在进行性能分析时,您如何确保CPU的“Turbo Speed”功能对执行时间产生影响,即调制时钟频率? - codechimp
@Z boson FYI,代码的执行结果与预期不完全相同。由于A8被设置为全部为1,GCC会优化掉_mm256_mul_ps操作。需要将A8更改为其他值,例如5.0f。对于unroll1,缺少MULPS步骤不会影响总时间,但在unroll4中会有明显的影响。 - codechimp
@user4602856,有趣的观察。我大部分都是用MSVC完成的,可能没有优化掉1(我已经停止使用MSVC)。我可能会再次研究这个... - Z boson

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接