圆周率计算公式如下:
AVX相关指令介绍:https://software.intel.com/sites/landingpage/IntrinsicsGuide/
AVX指令中有计算反正切函数的指令,但本文的目的是为了对比验证AVX相关指令实现SIMD并行计算对于计算效率的提升。因此,我们选取定积分部分作为计算圆周率的衡量方式,即将计算圆周率视作大量的可迭代的加法运算。
本文相关代码:https://github.com/vistart/IntelIntrinsics
一、准备工作
1. 引用头文件
使用AVX相关协议需要引用头文件”immintrin.h”:
#include "immintrin.h"
2. 验证自己的处理器对AVX指令的支持程度
方法一:如果是英特尔处理器,且知道处理器的型号,可以到英特尔官方网站查询。
以英特尔酷睿i9-9900K为例,查询结果如下:

从图1中可以看出,英特尔酷睿i9-9900K支持AVX2指令集,但不支持AVX-512指令集。
再以英特尔酷睿i5-1135G7为例,查询结果如下:

从图2中可以看出,英特尔酷睿i5-1135G7支持AVX-512指令集。
方法二:如果不知道处理器的型号,或没有可查看处理器指令集支持情况的页面,可以下载CPU-Z后运行CPU-Z查询。CPU-Z会直接显示当前处理器支持的指令集情况。
以 AMD Ryzen 9 3900X 为例,指令集支持情况如下:

从图3中可以看出,AMD Ryzen 9 3900X支持AVX2指令集,不支持AVX-512指令集。
二、原始计算方式
原始计算方式是仅使用默认四则运算来计算圆周率。源代码如下:
double compute_pi_native(size_t dt) { double pi = 0.0; double const delta = 1.0 / dt; for (size_t i = 0; i < dt; i++) { double x = i * delta; pi += delta / (1 + x * x); } return pi * 4; }
三、使用AVX2相关指令计算
double compute_pi_avx256(size_t dt) { double pi = 0.0; double const delta = 1.0 / dt; __m256d ymm0, ymm1, ymm2, ymm3, ymm4; ymm0 = _mm256_set1_pd(1.0); ymm1 = _mm256_set1_pd(delta); ymm2 = _mm256_set_pd(delta * 3, delta * 2, delta, 0.0); ymm4 = _mm256_setzero_pd(); for (int i = 0; i < dt - 4; i += 4) { ymm3 = _mm256_set1_pd(i * delta); // x = i * delta ymm3 = _mm256_add_pd(ymm3, ymm2); // x1 = x, x2 = x + delta, x3 = x + 2 * delta, x4 = x + 3 * delta ymm3 = _mm256_mul_pd(ymm3, ymm3); // x * x ymm3 = _mm256_add_pd(ymm0, ymm3); // 1 + x * x ymm3 = _mm256_div_pd(ymm1, ymm3); // delta / (1 + x * x) ymm4 = _mm256_add_pd(ymm4, ymm3); // pi += delta / (1 + x * x) } double tmp[4]; _mm256_store_pd(tmp, ymm4); pi += tmp[0] + tmp[1] + tmp[2] + tmp[3]; return pi * 4.0; }
四、使用AVX-512相关指令计算
double compute_pi_avx512(size_t dt) { double pi = 0.0; double const delta = 1.0 / dt; __m512d ymm0, ymm1, ymm2, ymm3, ymm4; ymm0 = _mm512_set1_pd(1.0); ymm1 = _mm512_set1_pd(delta); ymm2 = _mm512_set_pd(delta * 7, delta * 6, delta * 5, delta * 4, delta * 3, delta * 2, delta, 0.0); ymm4 = _mm512_setzero_pd(); for (int i = 0; i < dt - 8; i += 8) { ymm3 = _mm512_set1_pd(i * delta); ymm3 = _mm512_add_pd(ymm3, ymm2); ymm3 = _mm512_mul_pd(ymm3, ymm3); ymm3 = _mm512_add_pd(ymm0, ymm3); ymm3 = _mm512_div_pd(ymm1, ymm3); ymm4 = _mm512_add_pd(ymm4, ymm3); } double tmp[8]; _mm512_store_pd(tmp, ymm4); pi += tmp[0] + tmp[1] + tmp[2] + tmp[3] + tmp[4] + tmp[5] + tmp[6] + tmp[7]; return pi * 4; }
五、实验
以下实验中,考量到每次运行耗时会有细微差异,每次Debug环境实验将运行十次,并计算耗时均值和方差;而 Release 模式会存在累进优化现象,即相同变量多次赋值只记录最后一次,因而 Release 下循环运行多次并分次计时会产生失真,故 Release 模式只通过观察方式记录十次耗时的均值和方差。
1. 实验环境一
编译环境: Visual Studio 2019,C++ 17
处理器:英特尔酷睿i5-1135G7,内存:16GB 4266MHz 四通道
在 Debug, x64 模式下,计算结果如下:
原始:3.14159265405554100070162348856684,30454毫秒
AVX2:3.14159265126151776570395668386482,15716毫秒
AVX-512:3.14159264753617062382318181334995,8362毫秒
在 Release, x64 模式下,计算结果如下:
原始:3.14159265405554100070162348856684,12325毫秒
AVX2:3.14159265126151776570395668386482,1045毫秒
AVX-512:3.14159264753617062382318181334995,1062毫秒
从以上结果可得出:Debug 模式下,AVX2指令的运行耗时比原始缩短了一半,而AVX-512又比AVX2缩短了一半。在 Release 模式下,AVX2指令的运行耗时相对原始缩短到十二分之一,而AVX-512与AVX2的运行耗时相差不大。反复测试若干次后,基本耗时和相互比例相差极小,可以认为此结果相对稳定。
2. 实验环境二
编译环境:Visual Studio 2019,C++ 17
处理器:英特尔酷睿 i9-9900K,内存:96GB,3000MHz 双通道
在 Debug, x64 模式下,计算结果如下:
原始:3.14159265405554100070162348856684,耗时均值18070.70毫秒,方差32.15毫秒
AVX2:3.14159265126151776570395668386482,耗时均值12475.30毫秒,方差35.33毫秒
在 Release, x64 模式下,计算结果如下:
原始:3.14159265405554100070162348856684,8650毫秒
AVX2:3.14159265126151776570395668386482,919毫秒
从以上结果可得出:Debug 模式下,AVX2指令的运行耗时约为原始的三分之二。在 Release 模式下,AVX2指令的运行还是低于原始的九分之一。因该处理器不支持AVX-512指令集,故此测试不涉及AVX-512指令集。
3. 实验环境三
编译环境:Visual Studio 2019,C++ 17
处理器:AMD Ryzen 3900X,内存:64GB,3000MHz 双通道
在 Debug, x64 模式下,计算结果如下:
原始:3.14159265405554100070162348856684,5770毫秒
AVX2:3.14159265126151776570395668386482,51343毫秒
在 Release, x64 模式下,计算结果如下:
原始:3.14159265405554100070162348856684,2616毫秒
AVX2:3.14159265126151776570395668386482,651毫秒
从以上结果可得出:Debug 模式下,AVX2指令的运行耗时接近原始的九倍。在 Release 模式下,AVX2指令的运行还是低于原始的四分之一。因该处理器不支持AVX-512指令集,故此测试不涉及AVX-512指令集。
六、总结
从实验结果中可看出,AVX指令集对并行计算效率提升十分显著。但在 Release 模式下 AVX2 和 AVX-512 未显示出明显差距,甚至 AVX-512 比 AVX2 还慢一点。不论是哪种方式,计算得出的结果有效位按十进制算,在小数点后7位。这个精度足以满足大多数计算需求,但更高精度计算仍待改进算法。接下来将深入分析其中原理,也欢迎广大网友提供建议。
您必须登录才能发表评论。