logo标签列表

Fast popcount on Intel Xeon Phi

cvectorizationxeon-phiintel-michammingweight
9
我正在Intel Xeon® Phi®上实现超快的popcount，因为它是各种生物信息学软件的性能热点。
我已经实现了五个代码片段。
#if defined(__MIC__)
#include <zmmintrin.h>
__attribute__((align(64))) static const uint32_t POPCOUNT_4bit[16] = {0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4};
__attribute__((align(64))) static const uint32_t MASK_4bit[16] = {0xF, 0xF, 0xF, 0xF, 0xF, 0xF, 0xF, 0xF, 0xF, 0xF, 0xF, 0xF, 0xF, 0xF, 0xF, 0xF};
inline uint64_t vpu_popcount1(uint64_t* buf, size_t n)  {
    register size_t result = 0;
    size_t i;
    register const __m512i popcnt = _mm512_load_epi32((void*)POPCOUNT_4bit);
    register const __m512i mask = _mm512_load_epi32((void*)MASK_4bit);
    register __m512i total;
    register __m512i shuf;

#pragma unroll(8)
    for (i = 0; i < n; i+=8) {
        shuf = _mm512_load_epi32(&buf[i]);
        _mm_prefetch((const char *)&buf[i+256], _MM_HINT_T1); // vprefetch1
        _mm_prefetch((const char *)&buf[i+64], _MM_HINT_T0); // vprefetch0
        total = _mm512_setzero_epi32();

        total = _mm512_add_epi32(_mm512_permutevar_epi32(_mm512_and_epi32(shuf, mask), popcnt), total);
        total = _mm512_add_epi32(_mm512_permutevar_epi32(_mm512_and_epi32(_mm512_srli_epi32(shuf, 4),  mask), popcnt), total);
        total = _mm512_add_epi32(_mm512_permutevar_epi32(_mm512_and_epi32(_mm512_srli_epi32(shuf, 8),  mask), popcnt), total);
        total = _mm512_add_epi32(_mm512_permutevar_epi32(_mm512_and_epi32(_mm512_srli_epi32(shuf, 12), mask), popcnt), total);
        total = _mm512_add_epi32(_mm512_permutevar_epi32(_mm512_and_epi32(_mm512_srli_epi32(shuf, 16), mask), popcnt), total);
        total = _mm512_add_epi32(_mm512_permutevar_epi32(_mm512_and_epi32(_mm512_srli_epi32(shuf, 20), mask), popcnt), total);
        total = _mm512_add_epi32(_mm512_permutevar_epi32(_mm512_and_epi32(_mm512_srli_epi32(shuf, 24), mask), popcnt), total);
        total = _mm512_add_epi32(_mm512_permutevar_epi32(_mm512_and_epi32(_mm512_srli_epi32(shuf, 28), mask), popcnt), total);

        /* Reduce add, which is analogous to SSSE3's PSADBW instruction,
           is not implementated as a single instruction in VPUv1, thus
           emulated by multiple instructions*/
        result += _mm512_reduce_add_epi32(total);
    }

    return result;
}

__attribute__((align(64))) static const unsigned magic[] = {\
        0x55555555, 0x55555555, 0x55555555, 0x55555555,\
        0x55555555, 0x55555555, 0x55555555, 0x55555555,\
        0x55555555, 0x55555555, 0x55555555, 0x55555555,\
        0x55555555, 0x55555555, 0x55555555, 0x55555555,\
        0x33333333, 0x33333333, 0x33333333, 0x33333333,\
        0x33333333, 0x33333333, 0x33333333, 0x33333333,\
        0x33333333, 0x33333333, 0x33333333, 0x33333333,\
        0x33333333, 0x33333333, 0x33333333, 0x33333333,\
        0x0F0F0F0F, 0x0F0F0F0F, 0x0F0F0F0F, 0x0F0F0F0F,\
        0x0F0F0F0F, 0x0F0F0F0F, 0x0F0F0F0F, 0x0F0F0F0F,\
        0x0F0F0F0F, 0x0F0F0F0F, 0x0F0F0F0F, 0x0F0F0F0F,\
        0x0F0F0F0F, 0x0F0F0F0F, 0x0F0F0F0F, 0x0F0F0F0F,\
        0x00FF00FF, 0x00FF00FF, 0x00FF00FF, 0x00FF00FF,\
        0x00FF00FF, 0x00FF00FF, 0x00FF00FF, 0x00FF00FF,\
        0x00FF00FF, 0x00FF00FF, 0x00FF00FF, 0x00FF00FF,\
        0x00FF00FF, 0x00FF00FF, 0x00FF00FF, 0x00FF00FF,\
        0x0000FFFF, 0x0000FFFF, 0x0000FFFF, 0x0000FFFF,\
        0x0000FFFF, 0x0000FFFF, 0x0000FFFF, 0x0000FFFF,\
        0x0000FFFF, 0x0000FFFF, 0x0000FFFF, 0x0000FFFF,\
        0x0000FFFF, 0x0000FFFF, 0x0000FFFF, 0x0000FFFF,\
            0x000000FF, 0x000000FF, 0x000000FF, 0x000000FF,\
            0x000000FF, 0x000000FF, 0x000000FF, 0x000000FF,\
            0x000000FF, 0x000000FF, 0x000000FF, 0x000000FF,\
            0x000000FF, 0x000000FF, 0x000000FF, 0x000000FF
    };

inline uint64_t vpu_popcount2(uint64_t* buf, size_t n)  {
    register size_t result = 0;
    size_t i;

    register const __m512i B0 = _mm512_load_epi32((void*)(magic+0));
    register const __m512i B1 = _mm512_load_epi32((void*)(magic+16));
    register const __m512i B2 = _mm512_load_epi32((void*)(magic+32));
    register const __m512i B3 = _mm512_load_epi32((void*)(magic+48));
    register const __m512i B4 = _mm512_load_epi32((void*)(magic+64));
    register __m512i total;
    register __m512i shuf;

#pragma unroll(8)
    for (i = 0; i < n; i+=8) {
        shuf = _mm512_load_epi32(&buf[i]);
        _mm_prefetch((const char *)&buf[i+512], _MM_HINT_T1); // vprefetch1
        _mm_prefetch((const char *)&buf[i+64], _MM_HINT_T0); // vprefetch0
        total = _mm512_sub_epi32(shuf, _mm512_and_epi32(B0, _mm512_srli_epi32(shuf,1)));
        total = _mm512_add_epi32(_mm512_and_epi32(B1, total), _mm512_and_epi32(B1,_mm512_srli_epi32(total,2)));
        total = _mm512_and_epi32(B2, _mm512_add_epi32(total, _mm512_srli_epi32(total,4)));
        total = _mm512_and_epi32(B3, _mm512_add_epi32(total, _mm512_srli_epi32(total,8)));
        total = _mm512_and_epi32(B4, _mm512_add_epi32(total, _mm512_srli_epi32(total,16)));

        /* Reduce add, which is analogous to SSSE3's PSADBW instruction,
           is not implementated as a single instruction in VPUv1, thus
           emulated by multiple instructions*/
        result += _mm512_reduce_add_epi32(total);
    }

    return result;
}

inline uint64_t vpu_popcount3(uint64_t* buf, size_t n)  {
    register size_t result = 0;
    size_t i;

    register const __m512i B0 = _mm512_load_epi32((void*)(magic+0));
    register const __m512i B1 = _mm512_load_epi32((void*)(magic+16));
    register const __m512i B2 = _mm512_load_epi32((void*)(magic+32));
    register const __m512i B3 = _mm512_load_epi32((void*)(magic+48));
    register const __m512i B4 = _mm512_load_epi32((void*)(magic+64));
    register __m512i total;
    register __m512i shuf;

#pragma unroll(4)
    for (i = 0; i < n; i+=16) {
        shuf = _mm512_load_epi32(&buf[i]);
        result += _mm_countbits_64(buf[i+8]);
        _mm_prefetch((const char *)&buf[i+512], _MM_HINT_T1); // vprefetch1
        _mm_prefetch((const char *)&buf[i+576], _MM_HINT_T1); // vprefetch1
        result += _mm_countbits_64(buf[i+9]);
        _mm_prefetch((const char *)&buf[i+64], _MM_HINT_T0); // vprefetch0
        _mm_prefetch((const char *)&buf[i+128], _MM_HINT_T0); // vprefetch0
        total = _mm512_sub_epi32(shuf, _mm512_and_epi32(B0, _mm512_srli_epi32(shuf,1)));
        result += _mm_countbits_64(buf[i+10]);
        total = _mm512_add_epi32(_mm512_and_epi32(B1, total), _mm512_and_epi32(B1,_mm512_srli_epi32(total,2)));
        result += _mm_countbits_64(buf[i+11]);
        total = _mm512_and_epi32(B2, _mm512_add_epi32(total, _mm512_srli_epi32(total,4)));
        result += _mm_countbits_64(buf[i+12]);
        total = _mm512_and_epi32(B3, _mm512_add_epi32(total, _mm512_srli_epi32(total,8)));
        result += _mm_countbits_64(buf[i+13]);
        total = _mm512_and_epi32(B4, _mm512_add_epi32(total, _mm512_srli_epi32(total,16)));
        result += _mm_countbits_64(buf[i+14]);

        /* Reduce add, which is analogous to SSSE3's PSADBW instruction,
           is not implementated as a single instruction in VPUv1, thus
           emulated by multiple instructions*/
        result += _mm512_reduce_add_epi32(total);
        result += _mm_countbits_64(buf[i+15]);
    }

    return result;
}

/* Using VPU or SSE's machine intrinsic, CPUs not supporting SIMD 
 * will use compiler's implementation, the speed of which depends */
static inline size_t scalar_popcountu(unsigned *buf, size_t n) {
  register size_t cnt = 0;
  size_t i;
#pragma vector always
#pragma unroll(8)
  for (i = 0; i < n; i++) {
    cnt += _mm_countbits_32(buf[i]);
    _mm_prefetch((const char *)&buf[i+512], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[i+64], _MM_HINT_T0); // vprefetch0
  }
  return cnt;
}

static inline size_t scalar_popcountlu(uint64_t *buf, size_t n) {
  register size_t cnt = 0;
  size_t i;
#pragma vector always
#pragma unroll(8)
  for (i = 0; i < n; i++) {
    cnt += _mm_countbits_64(buf[i]);
    _mm_prefetch((const char *)&buf[i+512], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[i+64], _MM_HINT_T0); // vprefetch0
  }
  return cnt;
}
#endif

可以从https://www.dropbox.com/sh/b3sfqps19wa2oi4/iFQ9wQ1NTg下载带有OpenMP支持的代码。
使用Intel C/C++编译器XE 13编译该代码的命令如下：
icc -debug inline-debug-info -O3 -mmic -fno-alias -ansi-alias -opt-streaming-stores always -ipo popcnt-mmic.cpp -o popcnt-mmic -vec-report=2 -openmp

代码在协处理器上本地运行（61个内核），使用“122个线程”和线程亲和力为“平衡”，使用导出进行优化。
export OMP_NUM_THREADS=122;export KMP_AFFINITY=balanced

我正在使用Xeon Phi SE10p，B1步进，CentOS6.4 在28兆字节的垃圾数据（由rand()函数填充）上进行10000次迭代测试，性能如下：

Buffer allocated at: 0x7f456b000000
OpenMP scalar_popcountu       4310169 us; cnt = 28439328
OpenMP scalar_popcountlu      1421139 us; cnt = 28439328
OpenMP vpu_popcount           1489992 us; cnt = 28439328
OpenMP vpu_popcount2          1109530 us; cnt = 28439328
OpenMP vpu_popcount3           951122 us; cnt = 28439328

"scalar_popcountu" 和 "scalar_popcountlu" 分别使用 "_mm_countbits_32" 和 "_mm_countbits_64" 内在函数，这些函数利用标量 "popcnt" 指令。设置 "#pragma vector always" 要求编译器将加载和求和向量化为每次 16 个无符号整数或 8 个无符号长整数，尽管 popcount 本身仍然是标量指令。
vpu_popcount1 的实现类似于 SSSE3 popcount 实现 http://wm.ite.pl/articles/sse-popcount.html。然而，1）Xeon Phi 不支持整数的打包字节操作（最小值为双字，即 32 位），2）它不实现 "Packed sum of absolute difference" 指令（例如 SSSE3 中的 _mm_sad_epu8），因此通过四组 "vpermf32x4"、"vpaddd" 和 "movslq" 的组合来执行归约加法。因此，该实现生成的指令比原始的 SSSE3 版本多得多。"
vpu_popcount2的实现类似于SSE2 popcount实现（可以参考“Hacker's Delight”）。该实现生成的指令比vpu_popcount1少，速度约快30％。然而，繁琐的“reduce add”仍然无法避免。
vpu_popcount3的实现非常适用于Xeon Phi。通过向量和标量操作的混合，它比vpu_popcount2快约15％（标量操作在向量操作中的插入是我实现中的闲暇之作，可以根据编译器生成的汇编代码重新排列标量操作，但据我所知，预期的改进有限）。该改进基于以下观察结果：1）Xeon Phi是按顺序调度的，2）每个时钟周期可以发出两个标量指令或“1个向量+ 1个标量”指令。我将展开次数从8减少到4以避免寄存器文件饱和。
每个函数中从内存到L2提前8个循环，从L2到L1提前1个循环的显式预取已将L1命中率从0.38提高到0.994。

展开循环可以使性能提高约15%。这是违反直觉的，因为Xeon Phi采用顺序调度。但是展开使得icc编译器能够尽可能多地进行编译时调度。

我们还有更多技巧来提高性能吗？

Brian Nickerson提供了两个更快的代码片段。

OpenMP vpu_popcount2          1110737 us; cnt = 28439328
OpenMP vpu_popcount3           951459 us; cnt = 28439328
OpenMP vpu_popcount3_r         815126 us; cnt = 28439328
OpenMP vpu_popcount5           746852 us; cnt = 28439328

vpu_popcount3_revised：

inline uint64_t vpu_popcount3_revised(uint64_t* buf, size_t n) {
  _mm_prefetch((const char *)&buf[0], _MM_HINT_T0); // vprefetch0
  _mm_prefetch((const char *)&buf[8], _MM_HINT_T0); // vprefetch0
  _mm_prefetch((const char *)&buf[16], _MM_HINT_T1); // vprefetch1
  _mm_prefetch((const char *)&buf[24], _MM_HINT_T1); // vprefetch1
  _mm_prefetch((const char *)&buf[32], _MM_HINT_T1); // vprefetch1
  _mm_prefetch((const char *)&buf[40], _MM_HINT_T1); // vprefetch1
  _mm_prefetch((const char *)&buf[48], _MM_HINT_T1); // vprefetch1
  _mm_prefetch((const char *)&buf[56], _MM_HINT_T1); // vprefetch1
  _mm_prefetch((const char *)&buf[64], _MM_HINT_T1); // vprefetch1
  _mm_prefetch((const char *)&buf[72], _MM_HINT_T1); // vprefetch1
  _mm_prefetch((const char *)&buf[80], _MM_HINT_T1); // vprefetch1
  _mm_prefetch((const char *)&buf[88], _MM_HINT_T1); // vprefetch1
  _mm_prefetch((const char *)&buf[96], _MM_HINT_T1); // vprefetch1
  _mm_prefetch((const char *)&buf[104], _MM_HINT_T1); // vprefetch1
  _mm_prefetch((const char *)&buf[112], _MM_HINT_T1); // vprefetch1
  _mm_prefetch((const char *)&buf[120], _MM_HINT_T1); // vprefetch1
  register size_t result;
  size_t i;

  register const __m512i B0 = _mm512_load_epi32((void*)(magic+0));
  register const __m512i B1 = _mm512_load_epi32((void*)(magic+16));
  register const __m512i B2 = _mm512_load_epi32((void*)(magic+32));
  register const __m512i B3 = _mm512_load_epi32((void*)(magic+48));
  register const __m512i B4 = _mm512_load_epi32((void*)(magic+64));
  register __m512i total0;
  register __m512i total1;
  register __m512i shuf0;
  register __m512i shuf1;
  register __m512i result0;
  register __m512i result1;

  result0 = _mm512_setzero_epi32();
  result1 = _mm512_setzero_epi32();

  for (i = 0; i < n; i+=16) {
      shuf0 = _mm512_load_epi32(&buf[i  ]);
      shuf1 = _mm512_load_epi32(&buf[i+8]);
      _mm_prefetch((const char *)&buf[i+128], _MM_HINT_T1); // vprefetch1
      _mm_prefetch((const char *)&buf[i+136], _MM_HINT_T1); // vprefetch1
      _mm_prefetch((const char *)&buf[i+16], _MM_HINT_T0); // vprefetch0
      _mm_prefetch((const char *)&buf[i+24], _MM_HINT_T0); // vprefetch0
      total0 = _mm512_sub_epi32(shuf0, _mm512_and_epi32(B0, _mm512_srli_epi32(shuf0,1)));
      total1 = _mm512_sub_epi32(shuf1, _mm512_and_epi32(B0, _mm512_srli_epi32(shuf1,1)));
      total0 = _mm512_add_epi32(_mm512_and_epi32(B1, total0), _mm512_and_epi32(B1,_mm512_srli_epi32(total0,2)));
      total1 = _mm512_add_epi32(_mm512_and_epi32(B1, total1), _mm512_and_epi32(B1,_mm512_srli_epi32(total1,2)));
      total0 = _mm512_and_epi32(B2, _mm512_add_epi32(total0, _mm512_srli_epi32(total0,4)));
      total1 = _mm512_and_epi32(B2, _mm512_add_epi32(total1, _mm512_srli_epi32(total1,4)));
      total0 = _mm512_and_epi32(B3, _mm512_add_epi32(total0, _mm512_srli_epi32(total0,8)));
      total1 = _mm512_and_epi32(B3, _mm512_add_epi32(total1, _mm512_srli_epi32(total1,8)));
      total0 = _mm512_and_epi32(B4, _mm512_add_epi32(total0, _mm512_srli_epi32(total0,16)));
      total1 = _mm512_and_epi32(B4, _mm512_add_epi32(total1, _mm512_srli_epi32(total1,16)));
      result0 = _mm512_add_epi32(result0,total0);
      result1 = _mm512_add_epi32(result1,total1);

  }

  result0 = _mm512_add_epi32(result0,result1);
  result  = _mm512_reduce_add_epi32(result0);

  return result;
}

vpu_popcount5：

inline uint64_t vpu_popcount5(uint64_t* buf, size_t n)  {
    _mm_prefetch((const char *)&buf[0], _MM_HINT_T0); // vprefetch0
    _mm_prefetch((const char *)&buf[8], _MM_HINT_T0); // vprefetch0
    _mm_prefetch((const char *)&buf[16], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[24], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[32], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[40], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[48], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[56], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[64], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[72], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[80], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[88], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[96], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[104], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[112], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[120], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[128], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[136], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[144], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[152], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[160], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[168], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[176], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[184], _MM_HINT_T1); // vprefetch1
    register size_t result;
    size_t i;

    register const __m512i B0 = _mm512_load_epi32((void*)(magic+0));
    register const __m512i B1 = _mm512_load_epi32((void*)(magic+16));
    register const __m512i B2 = _mm512_load_epi32((void*)(magic+32));
    register const __m512i B3 = _mm512_load_epi32((void*)(magic+48));
    register const __m512i B4 = _mm512_load_epi32((void*)(magic+64));
    register const __m512i B6 = _mm512_load_epi32((void*)(magic+80));
    register __m512i total0;
    register __m512i total1;
    register __m512i total2;
    register __m512i total3;
    register __m512i shuf0;
    register __m512i shuf1;
    register __m512i shuf2;
    register __m512i shuf3;
    register __m512i result0;
    register __m512i result1;

    result0 = _mm512_setzero_epi32();
    result1 = _mm512_setzero_epi32();

    for (i = 0; i < n; i+=32) {
            shuf0 = _mm512_load_epi32(&buf[i   ]);
            shuf1 = _mm512_load_epi32(&buf[i+ 8]);
            shuf2 = _mm512_load_epi32(&buf[i+16]);
            shuf3 = _mm512_load_epi32(&buf[i+24]);
            _mm_prefetch((const char *)&buf[i+192], _MM_HINT_T1); // vprefetch1
            _mm_prefetch((const char *)&buf[i+200], _MM_HINT_T1); // vprefetch1
            _mm_prefetch((const char *)&buf[i+208], _MM_HINT_T1); // vprefetch1
            _mm_prefetch((const char *)&buf[i+216], _MM_HINT_T1); // vprefetch1
            _mm_prefetch((const char *)&buf[i+32], _MM_HINT_T0); // vprefetch0
            _mm_prefetch((const char *)&buf[i+40], _MM_HINT_T0); // vprefetch0
            _mm_prefetch((const char *)&buf[i+48], _MM_HINT_T0); // vprefetch0
            _mm_prefetch((const char *)&buf[i+56], _MM_HINT_T0); // vprefetch0
            total0 = _mm512_sub_epi32(shuf0, _mm512_and_epi32(B0, _mm512_srli_epi32(shuf0,1)));                        //  max value in nn is 10
            total1 = _mm512_sub_epi32(shuf1, _mm512_and_epi32(B0, _mm512_srli_epi32(shuf1,1)));
            total2 = _mm512_sub_epi32(shuf2, _mm512_and_epi32(B0, _mm512_srli_epi32(shuf2,1)));
            total3 = _mm512_sub_epi32(shuf3, _mm512_and_epi32(B0, _mm512_srli_epi32(shuf3,1)));
            total0 = _mm512_add_epi32(_mm512_and_epi32(B1, total0), _mm512_and_epi32(B1,_mm512_srli_epi32(total0,2))); //  max value in nnnn is 0100
            total1 = _mm512_add_epi32(_mm512_and_epi32(B1, total1), _mm512_and_epi32(B1,_mm512_srli_epi32(total1,2)));
            total2 = _mm512_add_epi32(_mm512_and_epi32(B1, total2), _mm512_and_epi32(B1,_mm512_srli_epi32(total2,2)));
            total3 = _mm512_add_epi32(_mm512_and_epi32(B1, total3), _mm512_and_epi32(B1,_mm512_srli_epi32(total3,2)));
            total0 = _mm512_and_epi32(B2, _mm512_add_epi32(total0, _mm512_srli_epi32(total0,4)));                      //  max value in 0000nnnn is 00001000
            total1 = _mm512_and_epi32(B2, _mm512_add_epi32(total1, _mm512_srli_epi32(total1,4)));
            total2 = _mm512_and_epi32(B2, _mm512_add_epi32(total2, _mm512_srli_epi32(total2,4)));
            total3 = _mm512_and_epi32(B2, _mm512_add_epi32(total3, _mm512_srli_epi32(total3,4)));
            total0 = _mm512_add_epi32(total0, total1);                                                                 //  max value in 000nnnnn is 00010000
            total1 = _mm512_add_epi32(total2, total3);
            total0 = _mm512_add_epi32(total0, _mm512_srli_epi32(total0,8));                                            //  max value in xxxxxxxx00nnnnnn is 00100000
            total1 = _mm512_add_epi32(total1, _mm512_srli_epi32(total1,8));
            total0 = _mm512_and_epi32(B6, _mm512_add_epi32(total0, _mm512_srli_epi32(total0,16)));                     //  max value in each element is 01000000, i.e. 64
            total1 = _mm512_and_epi32(B6, _mm512_add_epi32(total1, _mm512_srli_epi32(total1,16)));
            result0 = _mm512_add_epi32(result0,total0);
            result1 = _mm512_add_epi32(result1,total1);
    }

    result0 = _mm512_add_epi32(result0,result1);
    result  = _mm512_reduce_add_epi32(result0);

    return result;
}
2
我认为这个问题应该放在Code Review上，因为你的代码似乎已经可以工作了。 - Morwenn
1
非常有趣的问题 - 我希望您能在社区中保持活跃。我想知道您用来测量缓存命中率的技术是什么。干杯！ - BlueStrat
@BlueStrat，我使用英特尔的vTune来获取缓存命中率。 - Aquaskyline
2个回答
2
自从昨天发布以来，我已经能够在我的卡片上运行您的代码和我的建议。由于硬件的步进和编译器版本可能存在关联，因此我的计时结果与您的不完全相同。但是总体趋势保持不变，我的建议似乎实现了约15%的性能提升。
通过稍微调整以下代码，我获得了额外的小幅度提升，大约在5%至10%之间。请注意，在以下代码片段中，B6的每个元素都设置为0x000000FF。此时，我认为该算法可能已经接近从GDDR到L2缓存可持续提供的最大带宽。
（附注：这种说法的证明是，如果我用一个for循环将popcount5函数的主体重复十次 - 请注意，这是对“chunk_size”输入数据的十次快速重复，所以其中九次它将非常热门在L2中 - 测试的总时间仅增加了约五倍，而不是十倍。我提出这个问题是因为我认为您的目标是调整位计数逻辑的速度，但是您希望部署它的应用程序实际上具有更小和/或更热的工作集。如果是这样，DRAM->L2带宽引入的节流正在模糊图片。但是请注意，将测试输入的大小降低，以使其保持在L2中变得更加热门似乎会导致其他开销 - 可能是openmp开销 - 相对更显着。）
inline uint64_t vpu_popcount5(uint64_t* buf, size_t n)  {
    _mm_prefetch((const char *)&buf[0], _MM_HINT_T0); // vprefetch0
    _mm_prefetch((const char *)&buf[8], _MM_HINT_T0); // vprefetch0
    _mm_prefetch((const char *)&buf[16], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[24], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[32], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[40], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[48], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[56], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[64], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[72], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[80], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[88], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[96], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[104], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[112], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[120], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[128], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[136], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[144], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[152], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[160], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[168], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[176], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[184], _MM_HINT_T1); // vprefetch1
    register size_t result;
    size_t i;

    register const __m512i B0 = _mm512_load_epi32((void*)(magic+0));
    register const __m512i B1 = _mm512_load_epi32((void*)(magic+16));
    register const __m512i B2 = _mm512_load_epi32((void*)(magic+32));
    register const __m512i B6 = _mm512_load_epi32((void*)(magic+80));
    register __m512i total0;
    register __m512i total1;
    register __m512i total2;
    register __m512i total3;
    register __m512i shuf0;
    register __m512i shuf1;
    register __m512i shuf2;
    register __m512i shuf3;
    register __m512i result0;
    register __m512i result1;

    result0 = _mm512_setzero_epi32();
    result1 = _mm512_setzero_epi32();

    for (i = 0; i < n; i+=32) {
        shuf0 = _mm512_load_epi32(&buf[i   ]);
        shuf1 = _mm512_load_epi32(&buf[i+ 8]);
        shuf2 = _mm512_load_epi32(&buf[i+16]);
        shuf3 = _mm512_load_epi32(&buf[i+24]);
        _mm_prefetch((const char *)&buf[i+192], _MM_HINT_T1); // vprefetch1
        _mm_prefetch((const char *)&buf[i+200], _MM_HINT_T1); // vprefetch1
        _mm_prefetch((const char *)&buf[i+208], _MM_HINT_T1); // vprefetch1
        _mm_prefetch((const char *)&buf[i+216], _MM_HINT_T1); // vprefetch1
        _mm_prefetch((const char *)&buf[i+32], _MM_HINT_T0); // vprefetch0
        _mm_prefetch((const char *)&buf[i+40], _MM_HINT_T0); // vprefetch0
        _mm_prefetch((const char *)&buf[i+48], _MM_HINT_T0); // vprefetch0
        _mm_prefetch((const char *)&buf[i+56], _MM_HINT_T0); // vprefetch0
        total0 = _mm512_sub_epi32(shuf0, _mm512_and_epi32(B0, _mm512_srli_epi32(shuf0,1)));                        //  max value in nn is 10
        total1 = _mm512_sub_epi32(shuf1, _mm512_and_epi32(B0, _mm512_srli_epi32(shuf1,1)));
        total2 = _mm512_sub_epi32(shuf2, _mm512_and_epi32(B0, _mm512_srli_epi32(shuf2,1)));
        total3 = _mm512_sub_epi32(shuf3, _mm512_and_epi32(B0, _mm512_srli_epi32(shuf3,1)));
        total0 = _mm512_add_epi32(_mm512_and_epi32(B1, total0), _mm512_and_epi32(B1,_mm512_srli_epi32(total0,2))); //  max value in nnnn is 0100
        total1 = _mm512_add_epi32(_mm512_and_epi32(B1, total1), _mm512_and_epi32(B1,_mm512_srli_epi32(total1,2)));
        total2 = _mm512_add_epi32(_mm512_and_epi32(B1, total2), _mm512_and_epi32(B1,_mm512_srli_epi32(total2,2)));
        total3 = _mm512_add_epi32(_mm512_and_epi32(B1, total3), _mm512_and_epi32(B1,_mm512_srli_epi32(total3,2)));
        total0 = _mm512_and_epi32(B2, _mm512_add_epi32(total0, _mm512_srli_epi32(total0,4)));                      //  max value in 0000nnnn is 00001000
        total1 = _mm512_and_epi32(B2, _mm512_add_epi32(total1, _mm512_srli_epi32(total1,4)));
        total2 = _mm512_and_epi32(B2, _mm512_add_epi32(total2, _mm512_srli_epi32(total2,4)));
        total3 = _mm512_and_epi32(B2, _mm512_add_epi32(total3, _mm512_srli_epi32(total3,4)));
        total0 = _mm512_add_epi32(total0, total1);                                                                 //  max value in 000nnnnn is 00010000
        total1 = _mm512_add_epi32(total2, total3);
        total0 = _mm512_add_epi32(total0, _mm512_srli_epi32(total0,8));                                            //  max value in xxxxxxxx00nnnnnn is 00100000
        total1 = _mm512_add_epi32(total1, _mm512_srli_epi32(total1,8));
        total0 = _mm512_and_epi32(B6, _mm512_add_epi32(total0, _mm512_srli_epi32(total0,16)));                     //  max value in each element is 01000000, i.e. 64
        total1 = _mm512_and_epi32(B6, _mm512_add_epi32(total1, _mm512_srli_epi32(total1,16)));
        result0 = _mm512_add_epi32(result0,total0);
        result1 = _mm512_add_epi32(result1,total1);

        /* Reduce add, which is analogous to SSSE3's PSADBW instruction,
           is not implementated as a single instruction in VPUv1, thus
           emulated by multiple instructions*/
    }

    result0 = _mm512_add_epi32(result0,result1);
    result  = _mm512_reduce_add_epi32(result0);

    return result;
}
代码很有启发性，但结果不正确。我已经将B6的使用更改回B3和B4，然后再次获得了正确的结果。 B6： OpenMP vpu_popcount3 918206微秒; cnt = 28439328 OpenMP vpu_popcount3_r 782216微秒; cnt = 28439328 OpenMP vpu_popcount5 751597微秒; cnt = 3582282 B3，B4： OpenMP vpu_popcount3 951459微秒; cnt = 28439328 OpenMP vpu_popcount3_r 815126微秒; cnt = 28439328 OpenMP vpu_popcount5 746852微秒; cnt = 28439328 - Aquaskyline
我正在使用B1步进，SE10p。 - Aquaskyline
我仔细检查了一遍，我得到了正确的答案。你是否在“magic”数组的末尾添加了16个0x000000FF的副本？（我在代码片段之前的散文中提到过这一点，但我可能应该更直接地说明。抱歉。请注意，为了与其他Bn寄存器常量的命名一致，我应该将其称为B5，但这并不重要。） - Brian Nickerson
我在度假前进行了最后一个实验：我黑进了popcount_openmp函数中对popcount函数的调用，使其总是传递buf[0]的地址而不是buf[i]的地址。这样做的原因是为了保证数据在L2中非常热，以便我能够确定我对L2带宽的担忧是否成立。显然并没有；时间并没有很大的差异。另一个实验的结果（重复函数体10次只需要5倍的时间）可能是由于相对较少的openmp开销。 - Brian Nickerson
已更新代码，加入了新的“神奇”数组，现在可以正常工作了。 - Aquaskyline
1
请您尝试以下变体，并回报是否对性能有所改善？我正在解决我认为您编码中不太优化的几个问题：
- 我认为您的预取距离不太正确。在索引实际上是基于uint64时，看起来您可能在考虑字节偏移距离。 - 我认为没有必要在每次循环迭代中执行缩减操作。您可以在16个SIMD元素中进行位数部分累加，然后在循环外部执行单个缩减操作。 - 我认为执行标量端popcount指令并不如真正获得最佳VPU调度优势。专注于优秀的VPU调度非常重要。我还认为标量popcount指令实际上并不与矢量操作配对；即我认为它仅在U管道中支持。 
inline uint64_t vpu_popcount3_revised(uint64_t* buf, size_t n) {
    _mm_prefetch((const char *)&buf[0], _MM_HINT_T0); // vprefetch0
    _mm_prefetch((const char *)&buf[8], _MM_HINT_T0); // vprefetch0
    _mm_prefetch((const char *)&buf[16], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[24], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[32], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[40], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[48], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[56], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[64], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[72], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[80], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[88], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[96], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[104], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[112], _MM_HINT_T1); // vprefetch1
    _mm_prefetch((const char *)&buf[120], _MM_HINT_T1); // vprefetch1
    register size_t result;
    size_t i;

    register const __m512i B0 = _mm512_load_epi32((void*)(magic+0));
    register const __m512i B1 = _mm512_load_epi32((void*)(magic+16));
    register const __m512i B2 = _mm512_load_epi32((void*)(magic+32));
    register const __m512i B3 = _mm512_load_epi32((void*)(magic+48));
    register const __m512i B4 = _mm512_load_epi32((void*)(magic+64));
    register __m512i total0;
    register __m512i total1;
    register __m512i shuf0;
    register __m512i shuf1;
    register __m512i result0;
    register __m512i result1;

    result0 = _mm512_setzero_epi32();
    result1 = _mm512_setzero_epi32();

    for (i = 0; i < n; i+=16) {
        shuf0 = _mm512_load_epi32(&buf[i  ]);
        shuf1 = _mm512_load_epi32(&buf[i+8]);
        _mm_prefetch((const char *)&buf[i+128], _MM_HINT_T1); // vprefetch1
        _mm_prefetch((const char *)&buf[i+136], _MM_HINT_T1); // vprefetch1
        _mm_prefetch((const char *)&buf[i+16], _MM_HINT_T0); // vprefetch0
        _mm_prefetch((const char *)&buf[i+24], _MM_HINT_T0); // vprefetch0
        total0 = _mm512_sub_epi32(shuf0, _mm512_and_epi32(B0, _mm512_srli_epi32(shuf0,1)));
        total1 = _mm512_sub_epi32(shuf1, _mm512_and_epi32(B0, _mm512_srli_epi32(shuf1,1)));
        total0 = _mm512_add_epi32(_mm512_and_epi32(B1, total0), _mm512_and_epi32(B1,_mm512_srli_epi32(total0,2)));
        total1 = _mm512_add_epi32(_mm512_and_epi32(B1, total1), _mm512_and_epi32(B1,_mm512_srli_epi32(total1,2)));
        total0 = _mm512_and_epi32(B2, _mm512_add_epi32(total0, _mm512_srli_epi32(total0,4)));
        total1 = _mm512_and_epi32(B2, _mm512_add_epi32(total1, _mm512_srli_epi32(total1,4)));
        total0 = _mm512_and_epi32(B3, _mm512_add_epi32(total0, _mm512_srli_epi32(total0,8)));
        total1 = _mm512_and_epi32(B3, _mm512_add_epi32(total1, _mm512_srli_epi32(total1,8)));
        total0 = _mm512_and_epi32(B4, _mm512_add_epi32(total0, _mm512_srli_epi32(total0,16)));
        total1 = _mm512_and_epi32(B4, _mm512_add_epi32(total1, _mm512_srli_epi32(total1,16)));
        result0 = _mm512_add_epi32(result0,total0);
        result1 = _mm512_add_epi32(result1,total1);

    }

    /* Reduce add, which is analogous to SSSE3's PSADBW instruction,
       is not implementated as a single instruction in VPUv1, thus
       emulated by multiple instructions*/

    result0 = _mm512_add_epi32(result0,result1);
    result  = _mm512_reduce_add_epi32(result0);

    return result;
}
OpenMP ci_popcount 3871901微秒; cnt = 28439328
OpenMP ci_popcountll 1336664微秒; cnt = 28439328
OpenMP vpu_popcount1 1364152微秒; cnt = 28439328
OpenMP vpu_popcount2 1052478微秒; cnt = 28439328
OpenMP vpu_popcount3 918206微秒; cnt = 28439328
OpenMP vpu_popcount3_r 782216微秒; cnt = 28439328 - Aquaskyline
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接