AVX2代码慢,没有AVX2

Question

我一直在努力开始使用AVX2指令而运气不好(这个功能列表很有帮助).最后,我得到了我的第一个程序编译和做我想要的.我必须做的程序需要两个,u_char并且复合了两个.本质上,我使用它来解码存储在相机的u_char数组中的数据,但我认为这与此问题无关.

获得double两者的过程u_char是:

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

where 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.我按照这里的说明进行了检查,我的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;
}

完整的例子可以在这里找到.

Answer 1

在定时间隔内这么少量的工作很难准确测量. cols = 80只有20个__m256d向量.

我SKYLAKE微架构的系统上测试程序跳来跳去之间9.53674e-07 s,1.19209e-06 s并0 s为次,通常更快与AVX2版本.(我_mm_pause()在另一个核心上运行繁忙循环以最大速度挂接所有核心.它是一个桌面i7-6700k,因此所有核心共享相同的核心时钟频率.)

gettimeofday显然没有足够的精确度来衡量任何短暂的东西. struct timeval使用秒和微 -seconds,不纳秒. 但我确实看到AVL2版本在Skylake上更快,编译时使用g++ -O3 -march=native.我没有Haswell可以测试.我SKYLAKE微架构是使用硬件P状态电源管理,所以即使我不挂的时间提前了CPU的频率,这将上升到最大值非常快.Haswell没有这个功能,所以这是你的事情可能奇怪的另一个原因.

如果您想测量挂钟时间(而不是核心时钟周期),请std::chrono像普通人一样使用. 使用C++ 11进行便携式计时代码的正确方法.

预热效果将占主导地位,并且你将std::vector::resize()时间间隔包含在内部.两个不同的std::vector<double>对象必须分别分配内存,因此第二个可能需要从操作系统获取新页面并且需要更长时间.也许第一个能够从自由列表中获取内存,如果之前main(或其他内容cout <<)做了一些临时分配,然后缩小或释放它.

这里有很多可能性:首先,有些人报告说Haswell上的前几微秒看到256位向量指令运行较慢,就像在Skylake上测量的Agner Fog一样.

可能CPU决定在第二个定时间隔(AVX2一个)期间加速到最大涡轮增压.在i7-4700MQ(2.4GHz Haswell)上可能需要20k个时钟周期.(英特尔上的周期丢失？rdtsc和CPU_CLK_UNHALTED.REF_TSC之间的不一致).

也许在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(单个打包)的吞吐量是vsqrtpdHaswell的两倍(packed-double)吞吐量.(https://agner.org/optimize/).它还可以更容易地x * approx_rsqrt(x)用作快速近似sqrt(x),可能使用Newton-Raphson迭代从~12位精度上升到~24(几乎同样精确_mm256_sqrt_ps).请参阅快速矢量化rsqrt和SSE/AVX的倒数,具体取决于精度.(如果你在同一个循环中有足够的工作要做,而你没有在分频器吞吐量上遇到瓶颈,那么实际的sqrt指令就可以了.)

你可以SIMD sqrt与float和再转换为double,如果你真的需要你的输出格式是double为COMPAT你的代码的其余部分.

优化的东西其他比开方:

在Haswell上这可能不会更快,但如果其他线程没有使用SQRT/DIV,它可能更容易超线程.

它使用SIMD来加载和解包数据:a<<8 + b最好通过交换来自b和a生成16位整数的字节来实现_mm_unpacklo/hi_epi8.然后零扩展到32位整数,这样我们就可以使用SIMD int- > double转换.

这导致double每对__m128i数据的4个向量.这里使用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));
  }
}

这编译对Godbolt编译器Explorer中的一个相当不错的循环,用g++ -O3 -march=haswell.

要处理cols不是16的倍数,你需要另一个版本的循环,或填充或其他东西.

但是除了vsqrtpd根本没有任何指示之外没有任何帮助.

根据IACA的说法,Haswell瓶颈上的所有SIMD循环都在分配器单元上,每个循环28个循环vsqrtpd ymm,甚至是原始的,它会进行大量的标量工作.28个周期很长.

对于大输入,Skylake的速度应该是分频器吞吐量提高的两倍多.但是float仍然会加速~4倍,或更多vrsqrtps.