如何达到简单循环的 AVX 计算吞吐量

Question

最近我正在使用有限差分法研究计算电动力学的数值求解器。

求解器实现起来非常简单，但很难达到现代处理器的理论吞吐量，因为对加载的数据只有 1 次数学运算，例如：

    #pragma ivdep
    for(int ii=0;ii<Large_Number;ii++)
    { Z[ii]  = C1*Z[ii] + C2*D[ii];}

Large_Number 大约是 1,000,000，但不大于 10,000,000

我曾尝试手动展开循环并编写 AVX 代码，但未能使其更快：

int Vec_Size    = 8;
int Unroll_Num  = 6;
int remainder   = Large_Number%(Vec_Size*Unroll_Num);
int iter        = Large_Number/(Vec_Size*Unroll_Num);
int addr_incr   = Vec_Size*Unroll_Num;

__m256 AVX_Div1, AVX_Div2, AVX_Div3, AVX_Div4, AVX_Div5, AVX_Div6;
__m256 AVX_Z1,   AVX_Z2,   AVX_Z3,   AVX_Z4,   AVX_Z5,   AVX_Z6;

__m256 AVX_Zb = _mm256_set1_ps(Zb);
__m256 AVX_Za = _mm256_set1_ps(Za);
for(int it=0;it<iter;it++)
{
    int addr    = addr + addr_incr;
    AVX_Div1    = _mm256_loadu_ps(&Div1[addr]);     
    AVX_Z1      = _mm256_loadu_ps(&Z[addr]);
    AVX_Z1      = _mm256_add_ps(_mm256_mul_ps(AVX_Zb,AVX_Div1),_mm256_mul_ps(AVX_Za,AVX_Z1));
    _mm256_storeu_ps(&Z[addr],AVX_Z1);

    AVX_Div2    = _mm256_loadu_ps(&Div1[addr+8]);
    AVX_Z2      = _mm256_loadu_ps(&Z[addr+8]);
    AVX_Z2      = _mm256_add_ps(_mm256_mul_ps(AVX_Zb,AVX_Div2),_mm256_mul_ps(AVX_Za,AVX_Z2));
    _mm256_storeu_ps(&Z[addr+8],AVX_Z2);

    AVX_Div3    = _mm256_loadu_ps(&Div1[addr+16]);
    AVX_Z3      = _mm256_loadu_ps(&Z[addr+16]);
    AVX_Z3      = _mm256_add_ps(_mm256_mul_ps(AVX_Zb,AVX_Div3),_mm256_mul_ps(AVX_Za,AVX_Z3));
    _mm256_storeu_ps(&Z[addr+16],AVX_Z3);

    AVX_Div4    = _mm256_loadu_ps(&Div1[addr+24]);
    AVX_Z4      = _mm256_loadu_ps(&Z[addr+24]);
    AVX_Z4      = _mm256_add_ps(_mm256_mul_ps(AVX_Zb,AVX_Div4),_mm256_mul_ps(AVX_Za,AVX_Z4));
    _mm256_storeu_ps(&Z[addr+24],AVX_Z4);

    AVX_Div5    = _mm256_loadu_ps(&Div1[addr+32]);
    AVX_Z5      = _mm256_loadu_ps(&Z[addr+32]);
    AVX_Z5      = _mm256_add_ps(_mm256_mul_ps(AVX_Zb,AVX_Div5),_mm256_mul_ps(AVX_Za,AVX_Z5));
    _mm256_storeu_ps(&Z[addr+32],AVX_Z5);

    AVX_Div6    = _mm256_loadu_ps(&Div1[addr+40]);
    AVX_Z6      = _mm256_loadu_ps(&Z[addr+40]);
    AVX_Z6      = _mm256_add_ps(_mm256_mul_ps(AVX_Zb,AVX_Div6),_mm256_mul_ps(AVX_Za,AVX_Z6));
    _mm256_storeu_ps(&Z[addr+40],AVX_Z6);
}

上面的 AVX 循环实际上比 Inter 编译器生成的代码慢一点。

编译器生成的代码可以达到8G左右flops/s，大约是3GHz Ivybridge处理器单线程理论吞吐量的25%。我想知道是否有可能达到像这样的简单循环的吞吐量。

谢谢！

Answer 1

您确定 8 GFLOPS/s 大约是 3 GHz Ivybridge 处理器最大吞吐量的 25% 吗？让我们来计算一下。

每8个元素需要两次单精度AVX乘法和一次AVX加法。 Ivybridge 处理器每个周期只能执行一次 8 宽 AVX 加法和一次 8 宽 AVX 乘法。此外，由于加法依赖于两次乘法，因此需要 3 个周期来处理 8 个元素。由于加法可以与下一个乘法重叠，我们可以将其减少到每 8 个元素 2 个周期。对于十亿个元素，需要 2*10^9/8 = 10^9/4 个循环。考虑到 3 GHz 时钟，我们得到 10^9/4 * 10^-9/3 = 1/12 = 0.08 秒。所以最大理论吞吐量是12GLOPS/s，编译器生成的代码达到了66%，这很好。

还有一点，通过展开循环 8 次，可以有效地对其进行矢量化。我怀疑如果你展开这个特定的循环不止于此，尤其是超过 16 次，你是否会获得任何显着的加速。

Answer 2

我认为真正的瓶颈是每2个乘法和1个加法有2个加载和1个存储指令。也许计算是内存带宽限制。每个元素需要传输12字节数据，如果每秒处理2G元素（即6G flops）即24GB/s数据传输，达到ivy bridge的理论带宽。请问这个说法是否成立，确实没有解决这个问题的方法。

我之所以回答我自己的问题，是希望在我轻易放弃优化之前有人能指正我。这个简单的循环对于许多科学求解器来说极其重要，它是有限元和有限差分法的backbone。如果因为计算受内存带宽限制甚至不能满足一个处理器的需求，那为什么还要使用多核呢？高带宽 GPU 或 Xeon Phi 应该是更好的解决方案。

Answer 3

提高像您这样的代码的性能 "well explored" 并且仍然很受欢迎。查看点积（Z Boson 已经提供了完美的 link）或一些 (D)AXPY 优化讨论 (https://scicomp.stackexchange.com/questions/1932/are-daxpy-dcopy-dscal-overkills)

一般来说，探索和考虑应用的关键主题是：

由于 FMA 和更好的 load/store Haswell
预取。 "Streaming stores"、"non-temporal stores" 对于某些平台。
线程并行以达到最大持续带宽
展开（您已经完成；英特尔编译器也可以使用 #pragma unroll (X) 执行此操作）。 "streaming" 代码差别不大。
最终确定要针对哪些硬件平台优化代码

最后一个项目符号特别重要，因为对于 "streaming" 和整个内存绑定代码 - 了解更多关于目标内存系统的信息很重要；例如，对于现有的，尤其是未来的高端 HPC 服务器（例如代号为 Knights Landing 的第二代 Xeon Phi），您可能在带宽和计算之间有非常不同的 "roofline model" 平衡，甚至与优化情况下的技术不同对于普通台式机。

如何达到简单循环的 AVX 计算吞吐量

How to reach AVX computation throughput for a simple loop

c++

optimization

simd

intrinsics

avx