AVX2 代码比没有 AVX2 时慢
AVX2 code slower then without AVX2
我一直在尝试开始使用 AVX2 指令，但运气不佳（this 函数列表很有帮助）。最后，我得到了我的第一个程序编译并做我想做的事。我必须做的程序需要两个 u_char 并从中合成一个双倍。本质上，我用它来解码存储在相机 u_char 数组中的数据，但我认为与这个问题无关。
获取两个u_char的double的过程是：
double result = sqrt(double((msb<<8) + lsb)/64);
其中 msb 和 lsb 是具有最高有效位 (msb) 和较低有效位 (lsb) 的两个 u_char 变量要计算的 double。数据存储在表示行主矩阵的数组中，其中值编码列 i 的 msb 和 lsb 分别位于第二行和第三行。我已经使用和不使用 AVX2 对此进行了编码：
void getData(u_char* data, size_t cols, std::vector<double>& info)
{
info.resize(cols);
for (size_t i = 0; i < cols; i++)
{
info[i] = sqrt(double((data[cols + i] << 8) + data[2 * cols + i]) / 64.0);
;
}
}
void getDataAVX2(u_char* data, size_t cols, std::vector<double>& info)
{
__m256d dividend = _mm256_set_pd(1 / 64.0, 1 / 64.0, 1 / 64.0, 1 / 64.0);
info.resize(cols);
__m256d result;
for (size_t i = 0; i < cols / 4; i++)
{
__m256d divisor = _mm256_set_pd(double((data[4 * i + 3 + cols] << 8) + data[4 * i + 2 * cols + 3]),
double((data[4 * i + 2 + cols] << 8) + data[4 * i + 2 * cols + 2]),
double((data[4 * i + 1 + cols] << 8) + data[4 * i + 2 * cols + 1]),
double((data[4 * i + cols] << 8) + data[4 * i + 2 * cols]));
_mm256_storeu_pd(&info[0] + 4 * i, _mm256_sqrt_pd(_mm256_mul_pd(divisor, dividend)));
}
}
然而，令我惊讶的是，这段代码比正常代码慢？关于如何加快它的任何想法？
我正在使用 c++ (7.3.0) 使用以下选项 -std=c++17 -Wall -Wextra -O3 -fno-tree-vectorize -mavx2 进行编译。我已经按照 的说明进行了检查，我的 CPU（Intel(R) Core(TM) i7-4710HQ CPU @ 2.50GHz）支持 AVX2。
检查哪个更快是在使用时间。以下函数给我时间戳：
inline double timestamp()
{
struct timeval tp;
gettimeofday(&tp, nullptr);
return double(tp.tv_sec) + tp.tv_usec / 1000000.;
}
我在每个函数 getData 和 getDataAVX2 之前和之后获取时间戳，然后减去它们以获得每个函数的运行时间。总体main如下：
int main(int argc, char** argv)
{
u_char data[] = {
0xf, 0xf, 0xf, 0xf, 0xf, 0xf, 0xf, 0xf, 0xf, 0xf, 0x11, 0xf, 0xf, 0xf, 0xf, 0xf, 0x10, 0xf, 0xf,
0xf, 0xf, 0xe, 0x10, 0x10, 0xf, 0x10, 0xf, 0xf, 0x10, 0xf, 0xf, 0xf, 0xf, 0xf, 0xf, 0x10, 0x10, 0xf,
0x10, 0xf, 0xe, 0xf, 0xf, 0x10, 0xf, 0xf, 0x10, 0xf, 0xf, 0xf, 0xf, 0x10, 0xf, 0xf, 0xf, 0xf, 0xf,
0xf, 0xf, 0xf, 0x10, 0xf, 0xf, 0xf, 0x10, 0xf, 0xf, 0xf, 0xf, 0xe, 0xf, 0xf, 0xf, 0xf, 0xf, 0x10,
0x10, 0xf, 0xf, 0xf, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2,
0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2,
0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2,
0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2,
0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xd3, 0xd1, 0xca, 0xc6, 0xd2, 0xd2, 0xcc, 0xc8, 0xc2, 0xd0, 0xd0,
0xca, 0xc9, 0xcb, 0xc7, 0xc3, 0xc7, 0xca, 0xce, 0xca, 0xc9, 0xc2, 0xc8, 0xc2, 0xbe, 0xc2, 0xc0, 0xb8, 0xc4, 0xbd,
0xc5, 0xc9, 0xbc, 0xbf, 0xbc, 0xb5, 0xb6, 0xc1, 0xbe, 0xb7, 0xb9, 0xc8, 0xb9, 0xb2, 0xb2, 0xba, 0xb4, 0xb4, 0xb7,
0xad, 0xb2, 0xb6, 0xab, 0xb7, 0xaf, 0xa7, 0xa8, 0xa5, 0xaa, 0xb0, 0xa3, 0xae, 0xa9, 0xa0, 0xa6, 0xa5, 0xa8, 0x9f,
0xa0, 0x9e, 0x94, 0x9f, 0xa3, 0x9d, 0x9f, 0x9c, 0x9e, 0x99, 0x9a, 0x97, 0x4, 0x5, 0x4, 0x5, 0x4, 0x4, 0x5,
0x5, 0x5, 0x4, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x4, 0x4, 0x4, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5,
0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5,
0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x4, 0x4, 0x4, 0x5, 0x5, 0x5, 0x4, 0x4,
0x5, 0x5, 0x5, 0x5, 0x4, 0x5, 0x5, 0x4, 0x4, 0x6, 0x4, 0x4, 0x6, 0x5, 0x4, 0x5, 0xf0, 0xf0, 0xf0,
0xf0, 0xf0, 0xf0, 0xe0, 0xf0, 0xe0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0,
0xf0, 0xf0, 0xe0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0,
0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0,
0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0,
0xf0
};
size_t cols = 80;
// Normal
std::cout << "Computing with normal way" << std::endl;
std::vector<double> info;
double tstart_normal = timestamp();
getData(data, cols, info);
double time_normal = timestamp() - tstart_normal;
// AVX2
std::cout << "Computing with avx" << std::endl;
std::vector<double> info_avx2;
double tstart_avx2 = timestamp();
getDataAVX2(data, cols, info_avx2);
double time_avx2 = timestamp() - tstart_avx2;
// Display difference
std::cout << "Time normal: " << time_normal << " s" << std::endl;
std::cout << "Time AVX2: " << time_avx2 << " s" << std::endl;
std::cout << "Time improvement AVX2: " << time_normal / time_avx2 << std::endl;
// Write to file
std::ofstream file;
file.open("out.csv");
for (size_t i = 0; i < cols; i++)
{
file << info[size_t(i)] << "," << info_avx2[size_t(i)];
file << std::endl;
}
file.close();
// Exit
return 0;
}
可以找到完整的示例 here。
时间间隔内如此微小的工作量很难准确测量。 cols = 80 只有 20 __m256d 个向量。
您在我的 Skylake 系统上的测试程序有时会在 9.53674e-07 s、1.19209e-06 s 和 0 s 之间跳动，AVX2 版本通常更快。 （我在另一个核心上有一个 _mm_pause() 忙循环 运行ning 以最大速度固定所有核心。它是桌面 i7-6700k，因此所有核心共享相同的核心时钟频率。）
gettimeofday 显然远不够精确，无法测量任何短的东西。 struct timeval 使用秒和 微 秒，而不是纳秒。 但我确实一致地看到 AVX2 版本在 Skylake 上更快，使用 [=19= 编译].我没有要测试的 Haswell。我的 Skylake 使用硬件 P-state 电源管理，所以即使我没有提前固定 CPU 频率，它也会很快上升到最大值 。 Haswell 没有这个功能，所以这是你的东西可能很奇怪的另一个原因。
如果你想测量挂钟时间 (), use std::chrono like a normal person. .
热身效果将占主导地位，并且您将 std::vector::resize() 包含在定时间隔内 。两个不同的 std::vector<double> 对象必须分别分配内存，因此第二个对象可能需要从 OS 获取新页面并且需要更长的时间。也许第一个能够从空闲列表中获取内存，如果 main 之前的东西（或 cout << 中的东西）做了一些临时分配然后 sh运行k 或释放它。
这里有很多可能性：首先，有人报告说在 Haswell 上看到 256 位向量指令 运行 前几微秒变慢，like Agner Fog measured on Skylake。
可能 CPU 决定在 第二个时间间隔（AVX2 时间间隔）期间 加速到最大涡轮增压。这在 i7-4700MQ（2.4GHz Haswell）上可能需要 20k 个时钟周期。 ().
也许在 write 系统调用（来自 cout <<）之后，TLB 未命中或分支未命中对第二个函数的伤害更大？ （在您的内核中启用 Spectre + Meltdown 缓解措施后，您应该期望代码在从系统调用返回后立即变慢 运行。）
因为你没有使用 -ffast-math，GCC 不会把你的标量 sqrt 变成 rsqrtss 的近似值，特别是因为它是 double 而不是 float。否则就可以解释了。
查看时间如何与问题大小成比例，以确保您的微基准测试是合理的，除非您尝试测量瞬态/预热效果，重复工作很多次。如果它没有优化，只需在定时间隔内围绕函数调用重复循环（而不是尝试加起来时间来自多个间隔）。检查编译器生成的 asm，或者至少检查时间是否与重复计数呈线性关系。您可以将函数 __attribute__((noinline,noclone)) 作为一种阻止优化器跨重复循环迭代进行优化的方法。
除了预热效果之外，您的 SIMD 版本应该是 Haswell 上标量的大约 2 倍。
除法单元上的标量和 SIMD 版本瓶颈，即使在合并到 __m256d 之前输入的标量计算效率低下也是如此。 Haswell 的 FP divide/sqrt 硬件只有 128 位宽（因此 vsqrtpd ymm 被分成两个 128 位的一半）。但是标量只利用了可能吞吐量的一半。
float 会给您带来 4 倍的吞吐量提升：每个 SIMD 向量的元素数量是原来的两倍，并且 vsqrtps（packed-single）的吞吐量是 vsqrtpd（packed-双）在哈斯韦尔。 (https://agner.org/optimize/)。它还可以更容易地使用 x * approx_rsqrt(x) 作为 sqrt(x) 的快速近似值，可能使用 Newton-Raphson 迭代从 ~12 位精度提高到 ~24（几乎与 [=40 一样准确） =]).
参见 。 （如果您在同一个循环中有足够多的工作要做，而且您没有在分频器吞吐量上遇到瓶颈，那么实际的 sqrt 指令可能会很好。）
如果您真的需要输出格式，您可以将 SIMD sqrt 与 float 和 然后 转换为 double double 以与您的其余代码兼容 .
优化东西 其他 比 sqrt:
这在 Haswell 上可能不会更快，但如果其他线程不使用 SQRT/DIV.
，它可能对超线程更友好
它使用 SIMD 加载和解压缩数据：a<<8 + b 最好通过交错 b 和 a 的字节来完成16 位整数，_mm_unpacklo/hi_epi8。然后零扩展到 32 位整数，这样我们就可以使用 SIMD int->double 转换。
每对 __m128i 数据产生 4 个 double 向量。此处使用 256 位向量只会引入车道交叉问题，并且由于 _mm256_cvtepi32_pd(__m128i) 的工作原理，需要提取到 128。
我改为直接在输出中使用 _mm256_storeu_pd，而不是希望 gcc 优化一次一个元素的分配。
我还注意到编译器在每次存储后都会重新加载 &info[0]，因为它的别名分析无法证明 _mm256_storeu_pd 只是修改向量数据，而不是控制块。所以我将基地址分配给编译器确定没有指向自身的 double* 局部变量。
#include <immintrin.h>
#include <vector>
inline
__m256d cvt_scale_sqrt(__m128i vi){
__m256d vd = _mm256_cvtepi32_pd(vi);
vd = _mm256_mul_pd(vd, _mm256_set1_pd(1./64.));
return _mm256_sqrt_pd(vd);
}
// assumes cols is a multiple of 16
// SIMD for everything before the multiple/sqrt as well
// but probably no speedup because this and others just bottleneck on that.
void getDataAVX2_vector_unpack(const u_char*__restrict data, size_t cols, std::vector<double>& info_vec)
{
info_vec.resize(cols); // TODO: hoist this out of the timed region
double *info = &info_vec[0]; // our stores don't alias the vector control-block
// but gcc doesn't figure that out, so read the pointer into a local
for (size_t i = 0; i < cols / 4; i+=4)
{
// 128-bit vectors because packed int->double expands to 256-bit
__m128i a = _mm_loadu_si128((const __m128i*)&data[4 * i + cols]); // 16 elements
__m128i b = _mm_loadu_si128((const __m128i*)&data[4 * i + 2*cols]);
__m128i lo16 = _mm_unpacklo_epi8(b,a); // a<<8 | b packed 16-bit integers
__m128i hi16 = _mm_unpackhi_epi8(b,a);
__m128i lo_lo = _mm_unpacklo_epi16(lo16, _mm_setzero_si128());
__m128i lo_hi = _mm_unpackhi_epi16(lo16, _mm_setzero_si128());
__m128i hi_lo = _mm_unpacklo_epi16(hi16, _mm_setzero_si128());
__m128i hi_hi = _mm_unpackhi_epi16(hi16, _mm_setzero_si128());
_mm256_storeu_pd(&info[4*(i + 0)], cvt_scale_sqrt(lo_lo));
_mm256_storeu_pd(&info[4*(i + 1)], cvt_scale_sqrt(lo_hi));
_mm256_storeu_pd(&info[4*(i + 2)], cvt_scale_sqrt(hi_lo));
_mm256_storeu_pd(&info[4*(i + 3)], cvt_scale_sqrt(hi_hi));
}
}
这个compiles to a pretty nice loop on the Godbolt compiler explorer，和g++ -O3 -march=haswell。
要处理 cols 不是 16 的倍数，您将需要另一个版本的循环，或填充或其他东西。
但是除了 vsqrtpd 之外，使用更少的指令对解决瓶颈根本没有帮助。
According to IACA，除法器单元上 Haswell 瓶颈上的所有 SIMD 循环，每个 vsqrtpd ymm 28 个循环，甚至是你的原始循环，它做了大量的标量工作。 28个周期是长时间。
对于大输入，Skylake 的速度应该是原来的两倍多，因为它改进了除法器吞吐量。但是 float 仍然是 ~4 倍的加速，或者更多 vrsqrtps.
我一直在尝试开始使用 AVX2 指令，但运气不佳（this 函数列表很有帮助）。最后，我得到了我的第一个程序编译并做我想做的事。我必须做的程序需要两个 u_char 并从中合成一个双倍。本质上，我用它来解码存储在相机 u_char 数组中的数据，但我认为与这个问题无关。
获取两个u_char的double的过程是：
double result = sqrt(double((msb<<8) + lsb)/64);
其中 msb 和 lsb 是具有最高有效位 (msb) 和较低有效位 (lsb) 的两个 u_char 变量要计算的 double。数据存储在表示行主矩阵的数组中，其中值编码列 i 的 msb 和 lsb 分别位于第二行和第三行。我已经使用和不使用 AVX2 对此进行了编码：
void getData(u_char* data, size_t cols, std::vector<double>& info)
{
info.resize(cols);
for (size_t i = 0; i < cols; i++)
{
info[i] = sqrt(double((data[cols + i] << 8) + data[2 * cols + i]) / 64.0);
;
}
}
void getDataAVX2(u_char* data, size_t cols, std::vector<double>& info)
{
__m256d dividend = _mm256_set_pd(1 / 64.0, 1 / 64.0, 1 / 64.0, 1 / 64.0);
info.resize(cols);
__m256d result;
for (size_t i = 0; i < cols / 4; i++)
{
__m256d divisor = _mm256_set_pd(double((data[4 * i + 3 + cols] << 8) + data[4 * i + 2 * cols + 3]),
double((data[4 * i + 2 + cols] << 8) + data[4 * i + 2 * cols + 2]),
double((data[4 * i + 1 + cols] << 8) + data[4 * i + 2 * cols + 1]),
double((data[4 * i + cols] << 8) + data[4 * i + 2 * cols]));
_mm256_storeu_pd(&info[0] + 4 * i, _mm256_sqrt_pd(_mm256_mul_pd(divisor, dividend)));
}
}
然而，令我惊讶的是，这段代码比正常代码慢？关于如何加快它的任何想法？
我正在使用 c++ (7.3.0) 使用以下选项 -std=c++17 -Wall -Wextra -O3 -fno-tree-vectorize -mavx2 进行编译。我已经按照
检查哪个更快是在使用时间。以下函数给我时间戳：
inline double timestamp()
{
struct timeval tp;
gettimeofday(&tp, nullptr);
return double(tp.tv_sec) + tp.tv_usec / 1000000.;
}
我在每个函数 getData 和 getDataAVX2 之前和之后获取时间戳，然后减去它们以获得每个函数的运行时间。总体main如下：
int main(int argc, char** argv)
{
u_char data[] = {
0xf, 0xf, 0xf, 0xf, 0xf, 0xf, 0xf, 0xf, 0xf, 0xf, 0x11, 0xf, 0xf, 0xf, 0xf, 0xf, 0x10, 0xf, 0xf,
0xf, 0xf, 0xe, 0x10, 0x10, 0xf, 0x10, 0xf, 0xf, 0x10, 0xf, 0xf, 0xf, 0xf, 0xf, 0xf, 0x10, 0x10, 0xf,
0x10, 0xf, 0xe, 0xf, 0xf, 0x10, 0xf, 0xf, 0x10, 0xf, 0xf, 0xf, 0xf, 0x10, 0xf, 0xf, 0xf, 0xf, 0xf,
0xf, 0xf, 0xf, 0x10, 0xf, 0xf, 0xf, 0x10, 0xf, 0xf, 0xf, 0xf, 0xe, 0xf, 0xf, 0xf, 0xf, 0xf, 0x10,
0x10, 0xf, 0xf, 0xf, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2,
0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2,
0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2,
0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2,
0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xf2, 0xd3, 0xd1, 0xca, 0xc6, 0xd2, 0xd2, 0xcc, 0xc8, 0xc2, 0xd0, 0xd0,
0xca, 0xc9, 0xcb, 0xc7, 0xc3, 0xc7, 0xca, 0xce, 0xca, 0xc9, 0xc2, 0xc8, 0xc2, 0xbe, 0xc2, 0xc0, 0xb8, 0xc4, 0xbd,
0xc5, 0xc9, 0xbc, 0xbf, 0xbc, 0xb5, 0xb6, 0xc1, 0xbe, 0xb7, 0xb9, 0xc8, 0xb9, 0xb2, 0xb2, 0xba, 0xb4, 0xb4, 0xb7,
0xad, 0xb2, 0xb6, 0xab, 0xb7, 0xaf, 0xa7, 0xa8, 0xa5, 0xaa, 0xb0, 0xa3, 0xae, 0xa9, 0xa0, 0xa6, 0xa5, 0xa8, 0x9f,
0xa0, 0x9e, 0x94, 0x9f, 0xa3, 0x9d, 0x9f, 0x9c, 0x9e, 0x99, 0x9a, 0x97, 0x4, 0x5, 0x4, 0x5, 0x4, 0x4, 0x5,
0x5, 0x5, 0x4, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x4, 0x4, 0x4, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5,
0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5,
0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x5, 0x4, 0x4, 0x4, 0x5, 0x5, 0x5, 0x4, 0x4,
0x5, 0x5, 0x5, 0x5, 0x4, 0x5, 0x5, 0x4, 0x4, 0x6, 0x4, 0x4, 0x6, 0x5, 0x4, 0x5, 0xf0, 0xf0, 0xf0,
0xf0, 0xf0, 0xf0, 0xe0, 0xf0, 0xe0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0,
0xf0, 0xf0, 0xe0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0,
0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0,
0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0, 0xf0,
0xf0
};
size_t cols = 80;
// Normal
std::cout << "Computing with normal way" << std::endl;
std::vector<double> info;
double tstart_normal = timestamp();
getData(data, cols, info);
double time_normal = timestamp() - tstart_normal;
// AVX2
std::cout << "Computing with avx" << std::endl;
std::vector<double> info_avx2;
double tstart_avx2 = timestamp();
getDataAVX2(data, cols, info_avx2);
double time_avx2 = timestamp() - tstart_avx2;
// Display difference
std::cout << "Time normal: " << time_normal << " s" << std::endl;
std::cout << "Time AVX2: " << time_avx2 << " s" << std::endl;
std::cout << "Time improvement AVX2: " << time_normal / time_avx2 << std::endl;
// Write to file
std::ofstream file;
file.open("out.csv");
for (size_t i = 0; i < cols; i++)
{
file << info[size_t(i)] << "," << info_avx2[size_t(i)];
file << std::endl;
}
file.close();
// Exit
return 0;
}
可以找到完整的示例 here。
时间间隔内如此微小的工作量很难准确测量。 cols = 80 只有 20 __m256d 个向量。
您在我的 Skylake 系统上的测试程序有时会在 9.53674e-07 s、1.19209e-06 s 和 0 s 之间跳动，AVX2 版本通常更快。 （我在另一个核心上有一个 _mm_pause() 忙循环 运行ning 以最大速度固定所有核心。它是桌面 i7-6700k，因此所有核心共享相同的核心时钟频率。）
gettimeofday 显然远不够精确，无法测量任何短的东西。 struct timeval 使用秒和 微 秒，而不是纳秒。 但我确实一致地看到 AVX2 版本在 Skylake 上更快，使用 [=19= 编译].我没有要测试的 Haswell。我的 Skylake 使用硬件 P-state 电源管理，所以即使我没有提前固定 CPU 频率，它也会很快上升到最大值 。 Haswell 没有这个功能，所以这是你的东西可能很奇怪的另一个原因。
如果你想测量挂钟时间 (std::chrono like a normal person.
热身效果将占主导地位，并且您将 std::vector::resize() 包含在定时间隔内 。两个不同的 std::vector<double> 对象必须分别分配内存，因此第二个对象可能需要从 OS 获取新页面并且需要更长的时间。也许第一个能够从空闲列表中获取内存，如果 main 之前的东西（或 cout << 中的东西）做了一些临时分配然后 sh运行k 或释放它。
这里有很多可能性：首先，有人报告说在 Haswell 上看到 256 位向量指令 运行 前几微秒变慢，like Agner Fog measured on Skylake。
可能 CPU 决定在 第二个时间间隔（AVX2 时间间隔）期间 加速到最大涡轮增压。这在 i7-4700MQ（2.4GHz Haswell）上可能需要 20k 个时钟周期。 (
也许在 write 系统调用（来自 cout <<）之后，TLB 未命中或分支未命中对第二个函数的伤害更大？ （在您的内核中启用 Spectre + Meltdown 缓解措施后，您应该期望代码在从系统调用返回后立即变慢 运行。）
因为你没有使用 -ffast-math，GCC 不会把你的标量 sqrt 变成 rsqrtss 的近似值，特别是因为它是 double 而不是 float。否则就可以解释了。
查看时间如何与问题大小成比例，以确保您的微基准测试是合理的，除非您尝试测量瞬态/预热效果，重复工作很多次。如果它没有优化，只需在定时间隔内围绕函数调用重复循环（而不是尝试加起来时间来自多个间隔）。检查编译器生成的 asm，或者至少检查时间是否与重复计数呈线性关系。您可以将函数 __attribute__((noinline,noclone)) 作为一种阻止优化器跨重复循环迭代进行优化的方法。
除了预热效果之外，您的 SIMD 版本应该是 Haswell 上标量的大约 2 倍。
除法单元上的标量和 SIMD 版本瓶颈，即使在合并到 __m256d 之前输入的标量计算效率低下也是如此。 Haswell 的 FP divide/sqrt 硬件只有 128 位宽（因此 vsqrtpd ymm 被分成两个 128 位的一半）。但是标量只利用了可能吞吐量的一半。
float 会给您带来 4 倍的吞吐量提升：每个 SIMD 向量的元素数量是原来的两倍，并且 vsqrtps（packed-single）的吞吐量是 vsqrtpd（packed-双）在哈斯韦尔。 (https://agner.org/optimize/)。它还可以更容易地使用 x * approx_rsqrt(x) 作为 sqrt(x) 的快速近似值，可能使用 Newton-Raphson 迭代从 ~12 位精度提高到 ~24（几乎与 [=40 一样准确） =]).
参见
如果您真的需要输出格式，您可以将 SIMD sqrt 与 float 和 然后 转换为 double double 以与您的其余代码兼容 .
优化东西 其他 比 sqrt:
这在 Haswell 上可能不会更快，但如果其他线程不使用 SQRT/DIV.
，它可能对超线程更友好它使用 SIMD 加载和解压缩数据：a<<8 + b 最好通过交错 b 和 a 的字节来完成16 位整数，_mm_unpacklo/hi_epi8。然后零扩展到 32 位整数，这样我们就可以使用 SIMD int->double 转换。
每对 __m128i 数据产生 4 个 double 向量。此处使用 256 位向量只会引入车道交叉问题，并且由于 _mm256_cvtepi32_pd(__m128i) 的工作原理，需要提取到 128。
我改为直接在输出中使用 _mm256_storeu_pd，而不是希望 gcc 优化一次一个元素的分配。
我还注意到编译器在每次存储后都会重新加载 &info[0]，因为它的别名分析无法证明 _mm256_storeu_pd 只是修改向量数据，而不是控制块。所以我将基地址分配给编译器确定没有指向自身的 double* 局部变量。
#include <immintrin.h>
#include <vector>
inline
__m256d cvt_scale_sqrt(__m128i vi){
__m256d vd = _mm256_cvtepi32_pd(vi);
vd = _mm256_mul_pd(vd, _mm256_set1_pd(1./64.));
return _mm256_sqrt_pd(vd);
}
// assumes cols is a multiple of 16
// SIMD for everything before the multiple/sqrt as well
// but probably no speedup because this and others just bottleneck on that.
void getDataAVX2_vector_unpack(const u_char*__restrict data, size_t cols, std::vector<double>& info_vec)
{
info_vec.resize(cols); // TODO: hoist this out of the timed region
double *info = &info_vec[0]; // our stores don't alias the vector control-block
// but gcc doesn't figure that out, so read the pointer into a local
for (size_t i = 0; i < cols / 4; i+=4)
{
// 128-bit vectors because packed int->double expands to 256-bit
__m128i a = _mm_loadu_si128((const __m128i*)&data[4 * i + cols]); // 16 elements
__m128i b = _mm_loadu_si128((const __m128i*)&data[4 * i + 2*cols]);
__m128i lo16 = _mm_unpacklo_epi8(b,a); // a<<8 | b packed 16-bit integers
__m128i hi16 = _mm_unpackhi_epi8(b,a);
__m128i lo_lo = _mm_unpacklo_epi16(lo16, _mm_setzero_si128());
__m128i lo_hi = _mm_unpackhi_epi16(lo16, _mm_setzero_si128());
__m128i hi_lo = _mm_unpacklo_epi16(hi16, _mm_setzero_si128());
__m128i hi_hi = _mm_unpackhi_epi16(hi16, _mm_setzero_si128());
_mm256_storeu_pd(&info[4*(i + 0)], cvt_scale_sqrt(lo_lo));
_mm256_storeu_pd(&info[4*(i + 1)], cvt_scale_sqrt(lo_hi));
_mm256_storeu_pd(&info[4*(i + 2)], cvt_scale_sqrt(hi_lo));
_mm256_storeu_pd(&info[4*(i + 3)], cvt_scale_sqrt(hi_hi));
}
}
这个compiles to a pretty nice loop on the Godbolt compiler explorer，和g++ -O3 -march=haswell。
要处理 cols 不是 16 的倍数，您将需要另一个版本的循环，或填充或其他东西。
但是除了 vsqrtpd 之外，使用更少的指令对解决瓶颈根本没有帮助。
According to IACA，除法器单元上 Haswell 瓶颈上的所有 SIMD 循环，每个 vsqrtpd ymm 28 个循环，甚至是你的原始循环，它做了大量的标量工作。 28个周期是长时间。
对于大输入，Skylake 的速度应该是原来的两倍多，因为它改进了除法器吞吐量。但是 float 仍然是 ~4 倍的加速，或者更多 vrsqrtps.