c++ - AVX2 代码比没有 AVX2 慢

Question

我一直在尝试开始使用 AVX2 指令，但运气不佳(this 函数列表很有帮助)。最后，我编译了我的第一个程序并做我想做的事。我要做的程序需要两个u_char并把它加倍。本质上，我使用它来解码存储在来自相机的 u_char 数组中的数据，但我认为与此问题无关。

获取double的过程两者中的 u_char是:

double result = sqrt(double((msb<<8) + lsb)/64);

哪里msb和 lsb是两个 u_char具有最高有效位 ( msb ) 和较低有效位 ( lsb ) 的变量 double计算。数据存储在表示行主矩阵的数组中，其中 msb和 lsb值编码列 i分别在第二行和第三行。我在有和没有 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 .我已按照说明进行检查 here我的 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 中找到.

Answer 1

定时间隔内如此微量的工作量很难准确衡量。 cols = 80只有 20 __m256d vector 。

您在我的 Skylake 系统上的测试程序在 9.53674e-07 s 之间跳来跳去, 1.19209e-06 s和 0 s对于时代而言，使用 AVX2 版本通常会更快。 (我有一个 _mm_pause() 在另一个内核上运行的繁忙循环，以最大速度固定所有内核。它是台式机 i7-6700k，因此所有内核共享相同的内核时钟频率。)

gettimeofday显然远不足以测量任何这么短的东西。 struct timeval使用秒和微秒，而不是纳秒。 但我确实一直看到 AVX2 版本在 Skylake 上速度更快，编译通过 g++ -O3 -march=native .我没有要测试的 Haswell。我的 Skylake 使用硬件 P 状态电源管理，所以即使我没有提前调整 CPU 频率，它也会很快上升到最大值。 Haswell 没有这个功能，所以这是你的事情可能很奇怪的另一个原因。

如果要测量挂钟时间 ( instead of core clock cycles )，请使用 std::chrono像一个正常人。 Correct way of portably timing code using C++11 .

热身效果将占主导地位，您将包括 std::vector::resize()在定时间隔内 .两者不同std::vector<double>对象必须单独分配内存，因此第二个对象可能需要从操作系统获取新页面并花费更长的时间。也许第一个能够从空闲列表中获取内存，如果在 main 之前发生了什么| (或在 cout << 中的东西)做了一些临时分配，然后缩小或释放它。

这里有很多可能性:首先，有些人报告说在 Haswell 上看到 256 位 vector 指令在前几微秒内运行速度较慢，like Agner Fog measured on Skylake .

可能 CPU 决定在第二个定时间隔(AVX2 一个)内加速到最大涡轮增压。在 i7-4700MQ (2.4GHz Haswell) 上这可能需要 20k 个时钟周期。 ( Lost Cycles on Intel? An inconsistency between rdtsc and CPU_CLK_UNHALTED.REF_TSC )。

也许在 write 之后系统调用(来自 cout <<)TLB 未命中或分支未命中对第二个函数的伤害更大？ (在内核中启用 Spectre + Meltdown 缓解功能后，您应该期望代码在从系统调用返回后运行缓慢。)

由于您没有使用 -ffast-math , GCC 不会改变你的标量 sqrt成rsqrtss近似值，特别是因为它是 double不是 float .否则就可以解释了。

查看时间如何与问题大小成比例，以确保您的微基准是健全的，除非您尝试测量 transient /预热效果，否则请多次重复工作。 如果它没有优化掉，只需在定时间隔内的函数调用周围重复循环(而不是尝试将多个间隔的时间相加)。检查编译器生成的 asm，或者至少检查时间是否与重复计数成线性关系。您可能会创建函数 __attribute__((noinline,noclone))作为在重复循环迭代中优化优化器的一种方法。

在热身效果之外，您的 SIMD 版本应该是 Haswell 上标量的 2 倍左右 .

标量和 SIMD 版本都在除法单元上存在瓶颈，即使在合并到 __m256d 之前输入的标量计算效率低下. Haswell 的 FP 除法/sqrt 硬件只有 128 位宽(所以 vsqrtpd ymm 被分成两个 128 位的一半)。但是标量只利用了可能吞吐量的一半。
float将使您的吞吐量提高 4 倍:每个 SIMD vector 的元素数量增加两倍，并且 vsqrtps (packed-single) 的吞吐量是 vsqrtpd 的两倍(打包双)在 Haswell 上。 ( https://agner.org/optimize/ )。它还将使使用更容易x * approx_rsqrt(x)作为 sqrt(x) 的快速近似值，可能通过 Newton-Raphson 迭代从 ~12 位精度提高到 ~24(几乎与 _mm256_sqrt_ps 一样准确)。
见 Fast vectorized rsqrt and reciprocal with SSE/AVX depending on precision . (如果您有足够的工作在同一个循环中完成，并且您没有遇到分频器吞吐量的瓶颈，那么实际的 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转换。

这导致 double 的 4 个 vector 每对__m128i数据的。在这里使用 256 位 vector 只会引入车道交叉问题并且需要提取到 128，因为 _mm256_cvtepi32_pd(__m128i)作品。

我改用 _mm256_storeu_pd直接进入输出，而不是希望 gcc 会优化一次一个元素的分配。

我还注意到编译器正在重新加载 &info[0]在每家商店之后，因为它的别名分析无法证明_mm256_storeu_pd只是修改 vector 数据，而不是控制块。所以我将基地址分配给了 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 .