具有内在函数的条件代码的 C++ 向量化

Question

我尝试启用一个常用函数的矢量化来提高性能。

该算法应该执行以下操作并被调用 ~4.000.000 次！

Input:  double* cellvalue 
Output: int8*   Output (8 bit integer, c++ char) 

算法：

if (cellvalue > upper_threshold )
    *output = 1;
else if (cellvalue < lower_threshold)
    *output = -1;
else
    *output = 0;

我第一个并行计算 2 个双精度的矢量化方法如下所示：

__m128d lowerThresh = _mm_set1_pd(m_lowerThreshold);
__m128d upperThresh = _mm_set1_pd(m_upperThreshold);

__m128d vec = _mm_load_pd(cellvalue);
__m128d maskLower = _mm_cmplt_pd(vec, lowerThresh); // less than
__m128d maskUpper = _mm_cmpgt_pd(vec, upperThresh); // greater than

static const tInt8 negOne = -1;
static const tInt8 posOne =  1;
output[0] = (negOne & *((tInt8*)&maskLower.m128d_f64[0])) | (posOne & *((tInt8*)&maskUpper.m128d_f64[0]));
output[1] = (negOne & *((tInt8*)&maskLower.m128d_f64[1])) | (posOne & *((tInt8*)&maskUpper.m128d_f64[1]));

这对你有意义吗？它有效，但我认为创建输出的最后一部分非常复杂。有没有更快的方法来做到这一点？

我还尝试使用几乎相同的代码一次计算 8 个值。这会表现更好吗？指令的顺序有意义吗？

__m128d lowerThresh = _mm_set1_pd(m_lowerThreshold);
__m128d upperThresh = _mm_set1_pd(m_upperThreshold);

// load 4 times
__m128d vec0 = _mm_load_pd(cellValue);
__m128d vec1 = _mm_load_pd(cellValue + 2);
__m128d vec2 = _mm_load_pd(cellValue + 4);
__m128d vec3 = _mm_load_pd(cellValue + 6);
__m128d maskLower0 = _mm_cmplt_pd(vec0, lowerThresh); // less than
__m128d maskLower1 = _mm_cmplt_pd(vec1, lowerThresh); // less than
__m128d maskLower2 = _mm_cmplt_pd(vec2, lowerThresh); // less than
__m128d maskLower3 = _mm_cmplt_pd(vec3, lowerThresh); // less than
__m128d maskUpper0 = _mm_cmpgt_pd(vec0, upperThresh); // greater than
__m128d maskUpper1 = _mm_cmpgt_pd(vec1, upperThresh); // greater than
__m128d maskUpper2 = _mm_cmpgt_pd(vec2, upperThresh); // greater than
__m128d maskUpper3 = _mm_cmpgt_pd(vec3, upperThresh); // greater than

static const tInt8 negOne = -1;
static const tInt8 posOne =  1;
output[0] = (negOne & *((tInt8*)&maskLower0.m128d_f64[0])) | (posOne & *((tInt8*)&maskUpper0.m128d_f64[0]));
output[1] = (negOne & *((tInt8*)&maskLower0.m128d_f64[1])) | (posOne & *((tInt8*)&maskUpper0.m128d_f64[1]));
output[2] = (negOne & *((tInt8*)&maskLower1.m128d_f64[0])) | (posOne & *((tInt8*)&maskUpper1.m128d_f64[0]));
output[3] = (negOne & *((tInt8*)&maskLower1.m128d_f64[1])) | (posOne & *((tInt8*)&maskUpper1.m128d_f64[1]));
output[4] = (negOne & *((tInt8*)&maskLower2.m128d_f64[0])) | (posOne & *((tInt8*)&maskUpper2.m128d_f64[0]));
output[5] = (negOne & *((tInt8*)&maskLower2.m128d_f64[1])) | (posOne & *((tInt8*)&maskUpper2.m128d_f64[1]));
output[6] = (negOne & *((tInt8*)&maskLower3.m128d_f64[0])) | (posOne & *((tInt8*)&maskUpper3.m128d_f64[0]));
output[7] = (negOne & *((tInt8*)&maskLower3.m128d_f64[1])) | (posOne & *((tInt8*)&maskUpper3.m128d_f64[1]));

希望你能帮助我更好地理解矢量化的东西；）

Answer 1

_mm_cmplt_pd 和_mm_cmpgt_pd 产生的结果已经是0 或-1； and-1 什么都不做，and987654328@ 就等于否定它。因此，如果upper_threshold > lower_threshold（这样两个条件都不成立），你可以写*：

_mm_storeu_si128(output, _mm_sub_epi64(maskLower, maskUpper));

(*) 不清楚代码中的“int8”是什么；这不是 C++ 中的标准类型。它可能是一个 8 字节的 int，这是我在这里使用的行为。如果它是一个 8 位 int，则您需要打包一堆结果以存储在一起。

---

发问者澄清说，他们打算将int8 设为一个 8 位整数。在这种情况下，您可以执行以下操作以快速实现：

__m128i result = _mm_sub_epi64(maskLower, maskUpper)
output[0] = result.m128i_i64[0]; // .m128i_i64 is an oddball MSVC-ism, so
output[1] = result.m128i_i64[1]; // I'm not 100% sure about the syntax here.

但您可能还想尝试将八个结果向量打包在一起，并通过一次存储操作来存储它们。

Answer 2

如果您将代码更改为不分支，那么现代编译器将为您进行矢量化。

这是我运行的测试：

#include <stdint.h>                                                                             
#include <iostream>                                                                             
#include <random>                                                                               
#include <vector>                                                                               
#include <chrono>                                                                               

using Clock = std::chrono::steady_clock;                                                        
using std::chrono::milliseconds;                                                                

typedef double Scalar;                                                                          
typedef int8_t Integer;                                                                         

const Scalar kUpperThreshold = .5;                                                              
const Scalar kLowerThreshold = .2;                                                              

void compute_comparisons1(int n, const Scalar* xs, Integer* ys) {                               
#pragma simd                                                                                    
  for (int i=0; i<n; ++i) {                                                                     
    Scalar x   = xs[i];                                                                         
    ys[i] = (x > kUpperThreshold) - (x < kLowerThreshold);                                      
  }                                                                                             
}                                                                                               

void compute_comparisons2(int n, const Scalar* xs, Integer* ys) {                               
  for (int i=0; i<n; ++i) {                                                                     
    Scalar x   = xs[i];                                                                         
    Integer& y = ys[i];                                                                         
    if (x > kUpperThreshold)                                                                    
      y = 1;                                                                                    
    else if(x < kLowerThreshold)                                                                
      y = -1;                                                                                   
    else                                                                                        
      y = 0;                                                                                    
  }                                                                                             
}                                                                                               

const int N = 4000000;                                                                          

auto random_generator = std::mt19937{0};                                                        

int main() {                                                                                    
  std::vector<Scalar> xs(N);                                                                    
  std::vector<Integer> ys1(N);                                                                  
  std::vector<Integer> ys2(N);                                                                  

  std::uniform_real_distribution<Scalar> dist(0, 1);                                            
  for (int i=0; i<N; ++i)                                                                       
    xs[i] = dist(random_generator);                                                             


  auto time0 = Clock::now();                                                                    
  compute_comparisons1(N, xs.data(), ys1.data());                                               
  auto time1 = Clock::now();                                                                    
  compute_comparisons2(N, xs.data(), ys2.data());                                               
  auto time2 = Clock::now();                                                                    

  std::cout << "v1: " << std::chrono::duration_cast<milliseconds>(time1 - time0).count() << "\n";
  std::cout << "v2: " << std::chrono::duration_cast<milliseconds>(time2 - time1).count() << "\n";

  for (int i=0; i<N; ++i) {                                                                     
    if (ys1[i] != ys2[i]) {                                                                     
      std::cout << "Error!\n";                                                                  
      return -1;                                                                                
    }                                                                                           
  }                                                                                             
  return 0;                                                                                     
} 

如果您使用最新版本的 gcc（我使用 4.8.3）进行编译并使用标志“-O3 -std=c++11 -march=native -S”，您可以通过查看程序集来验证它将代码向量化。而且它的运行速度要快得多（在我的机器上是 3 毫秒对 16 毫秒。）

另外，我不确定您的要求是什么；但如果你能以较低的精度生活，那么使用 float 而不是 double 将进一步提高速度（在我的机器上 double 需要 1.8 倍）