如何实现简单循环的 AVX 计算吞吐量答案

【问题标题】：How to reach AVX computation throughput for a simple loop如何实现简单循环的 AVX 计算吞吐量
【发布时间】：2015-07-14 10:55:05
【问题描述】：

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

求解器实现起来非常简单，但很难达到现代处理器的理论吞吐量，因为对加载的数据只有 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%。我想知道是否有可能达到这样的简单循环的吞吐量。

谢谢！

【问题讨论】：

如果您的编译器支持，请在所有指针上尝试编译器特定的 restrict 关键字。如果不是，请在 C99 模式下编写和编译，并使用标准的restrict 关键字。此外，对齐所有指针，并告诉编译器它们是对齐的。另外，忘记余数，而是要求填充输入。
感谢您的 cmets，所有数组都使用 _mm_malloc(size,64) 与 64 字节对齐，并且使用 #pragma ivdep 而不是限制来向量化 for 循环。
你使用的是什么编译器和操作系统？
我看到一大堆未对齐的加载和存储（_mm256_loadu_ps、_mm256_storeu_ps）。是否可以不对齐数据以支持用对齐的版本（_mm256_load_ps、_mm256_store_ps）替换这些内在函数？

标签： c++ optimization simd intrinsics avx

【解决方案1】：

提高像您这样的代码的性能是“经过充分探索”并且仍然很受欢迎的领域。看看点积（Z Boson 已经提供了完美的链接）或一些 (D)AXPY 优化讨论 (https://scicomp.stackexchange.com/questions/1932/are-daxpy-dcopy-dscal-overkills)

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

由于 FMA 和 Haswell 上更好的加载/存储端口 u 架构，AVX2 优于 AVX
预取。某些平台的“流媒体商店”、“非临时商店”。
线程并行以达到最大持续带宽
展开（您已经完成；英特尔编译器也可以通过 #pragma unroll (X) 来完成）。 “流式传输”代码差别不大。
最后决定要优化代码的硬件平台是什么

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

【讨论】：

谢谢！我正在开发一个用于等离子体模拟的研究求解器，因为计算受到内存带宽的限制，只有 4 核 Xeon 的计算能力被利用了 8%。通过使用上述所有技巧来压缩最后一点带宽可能会将利用率提高到 10%（10 G flops/s 与 3GHz Ivy 的 100 Gf/s）。显然 CPU 不适合高带宽计算。
我想这就是英特尔决定不将 AVX512 包含在最新的 Skylake 处理器中的原因，只是没有足够的内存带宽来支持执行单元。
我不知道您针对哪些平台，但通常基于 Xeon 的平台具有更好的内存带宽绑定代码潜力。同样正如我所提到的，未来的 Xeons，尤其是 Xeon Phis 应该会显着改变平衡，有利于内存 bw 绑定的工作负载（即使存在 512 位宽的 SIMD）。
我没有看到任何关于“Skylakes 中不包括 AVX512”的公告。此外，还有足够多的受计算限制或延迟限制的工作负载，即使具有当前的内存带宽特性，它们（经过适当优化后）仍然可以从 AVX512 中受益；所以无论如何AVX512都不会白费。
AVX512 不包含在消费者 Skylakes 中，确切地说，它包含在 Xeon 版本中。

【解决方案2】：

您确定 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 秒。所以最大理论吞吐量是 12 GLOPS/s，编译器生成的代码达到了 66%，这很好。

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

【讨论】：

感谢您的回复。正如您所说，如果 8 个元素和 24 个触发器需要 2 个周期，那么实际上每个周期需要 12 个触发器，即 36G 触发器/秒。 25% 的估计或多或少是正确的
我认为真正的瓶颈是每 2 次乘法和 1 次加法都有 2 条加载和 1 条存储指令。也许计算是内存带宽有限的。每个元素需要传输 12 字节的数据，如果每秒处理 2G 个元素（即 6G flops）即 24GB/s 的数据传输，达到 ivy bridge 的理论带宽。我想知道这个论点是否成立，这个问题确实没有解决方案。
是的，你是对的。这相当于 36 GLOPS/s。
DRAM内存通道的数据速率是多少？
25.6GB/s 用于双通道 DDR3

【解决方案3】：

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

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

【讨论】：

您的循环受内存带宽限制。但是，这并不意味着您获得了峰值内存带宽。为此，您需要使用多个内核。这是使用多核的原因之一。见measuring-memory-bandwidth-from-the-dot-product-of-two-arrays。