c++ - 使用 Ivy Bridge 和 Haswell 循环展开以实现最大吞吐量

Question

我正在使用 AVX 一次计算八个点积。在我当前的代码中，我做了这样的事情(在展开之前):

Ivy 桥/沙桥

__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  ...

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

我猜这意味着不可能在 Ivy Bridge 和 Haswell 的 32 位模式下实现此任务的最大吞吐量，因为 32 位模式只有 8 个 AVX 寄存器？

编辑:关于赏金。我的主要问题仍然存在。我想获得上述 Ivy Bridge 或 Haswell 函数的最大吞吐量，n 可以是大于或等于 64 的任何值。我认为这只能使用展开来完成(Ivy 为八次桥牌和 Haswell 的 10 次)。如果你认为这可以用另一种方法来完成，那么让我们看看吧。在某种意义上，这是 How do I achieve the theoretical maximum of 4 FLOPs per cycle? 的变体。 .但不是只有乘法和加法，我正在寻找一个 256 位负载(或两个 128 位负载)、一个 AVX 乘法和一个 AVX 加法，每个时钟周期使用 Ivy Bridge 或两个 256 位负载和两个 FMA3 指令每个时钟周期。

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

我应该指出，编译器可能会给出误导性的结果 Difference in performance between MSVC and GCC for highly optimized matrix multplication code

我正在为 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);
}

Answer 1

对于 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 共同发出，因此它们不会影响展开因子。