使用Ivy Bridge和Haswell展开环路以实现最大吞吐量
Loop unrolling to achieve maximum throughput with Ivy Bridge and Haswell
我正在使用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);
}
Haswell
__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，我认为答案是沙桥和Haswell SSE2/AVX/AVX2的每个循环FLOPS。我需要把这个循环展开10次。
常春藤大桥我想是8号。这是我的逻辑。AVX加法具有3的延迟，乘法具有5的延迟。IvyBridge可以使用不同的端口同时进行一次AVX乘法和一次AVX加法。使用符号m表示乘法，a表示加法，x表示无运算，以及一个数字来表示部分和(例如，m5表示与第五个部分和相乘)，我可以写：
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位模式只有8个AVX寄存器
编辑：关于赏金。我的主要问题仍然存在。我想得到上面Ivy Bridge或Haswell函数的最大吞吐量，n可以是大于或等于64的任何值。我认为这只能通过展开来完成(Ivy Bridge为8次，Haswell为10次)。如果你认为这可以用另一种方法来完成，那么让我们看看。从某种意义上说，这是"如何实现每个周期4个FLOP的理论最大值？"的变体？。但是，我正在寻找一个256位加载(或两个128位加载)、一个AVX乘法和一个AVX加法，而不是只有乘法和加法，每个时钟周期使用Ivy Bridge或两个256位负载和两个FMA3指令。
我也想知道需要多少寄存器。对于常春藤大桥，我想是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);
}
对于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共同发布，因此它们不会影响展开因子
我在这里回答自己的问题只是为了添加信息。
我接着介绍了常春藤大桥的代码。当我第一次在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)
以下是i5-4250U在Linux(g++ -mavx -mfma -fopenmp -O3 main.cpp kernel_fma.cpp -o sum_fma)中使用带有GCC的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");
}
相关文章：
- 如果没有malloc，链表实现将失败
- 如何在c++中实现处理器调度模拟器
- 如何在c++中使用引用实现类似python的行为
- 实现无开销push_back的最佳方法是什么
- 使用简单类型列表实现的指数编译时间.为什么
- 如何在BST的这个简单递归实现中消除警告
- 实现一个在集合上迭代的模板函数
- 我应该实现右值推送功能吗?我应该使用std::move吗
- 如何正确实现和访问运算符的各种自定义枚举器
- C++Union/Struct位域的实现和可移植性
- 这个极客对极客的trie实现是否存在内存泄漏问题
- 在c++中实现LinkedList时，应出现未处理的错误
- 为左值和右值的包装器实现C++范围
- 使用模板进行堆栈实现; "name followed by :: must be a class or namespace"
- 使用GSoap实现ONVIF
- 在用于格式4的arm模拟器中实现功能时的一个问题
- 用于AVX的ln(x）的实现，m256
- 用常见虚拟函数实现的任意组合来实现派生类的正确方法是什么
- 使用Ivy Bridge和Haswell展开环路以实现最大吞吐量
- 在应用程序级别调整 TCP 数据包大小以实现最大吞吐量
