在C / C ++中使用AVX2的两个无符号字节向量的内积

编程问答 2022-06-07

问题描述

我想使用SSE / AVX2实现快速相关系数计算。操作数是两个unsigned char向量。该功能应大致等效于此：

float correlate_simple(const unsigned char* vec1,const unsigned char* vec2,size_t length)
{
    int sum1 = 0;
    int sum2 = 0;
    int sum11 = 0;
    int sum22 = 0;
    int sum12 = 0;

    for (size_t i = length; i > 0; --i,++vec1,++vec2) {
        sum1 += *vec1;
        sum2 += *vec2;
        sum11 += *vec1  * *vec1;
        sum22 += *vec2  * *vec2;
        sum12 += *vec1  * *vec2;
    }
    double mean1 = double(sum1) / double(length);
    double mean2 = double(sum2) / double(length);
    double mean11 = double(sum11) / double(length);
    double mean22 = double(sum22) / double(length);
    double mean12 = double(sum12) / double(length);

    double b = (mean11 - mean1 * mean1) * (mean22 - mean2 * mean2);
    if (b <= 0.0)
        return 0.0f;
    double a = (mean12 - mean1 * mean2);

    return float(a / sqrt(b));
}

参数length的范围是1到小于1000。

为此，我研究了如何实现两个无符号字节数组的内部乘积。但是，我无法提出一种不涉及将所有无符号的8位值转换为有符号的16位值的解决方案。

内部_mm256_maddubs_epi16(a,b)期望b是一个有符号字节。在这种情况下这将不是问题，因为从b中减去一些常数（此处为127）不会改变相关系数。不幸的是，我找不到能使我从无符号字节中减去127从而产生带符号字节的内在函数（不依赖于某些二进制补码的魔力）。

// vec: const unsigned char*
auto x = _mm256_load_si256((const __m256i*) vec);
auto v = _mm256_set1_epi8(127);

// wrong if vec[i] is less than 127:
auto x_centered = _mm256_sub_epi8 (x,v);

这里计算内部积（最后是相关系数）的最佳方法是什么？

附录：

下面是我目前对纯内部产品的实现。我决定转换为16位整数以避免溢出错误。 更新：从一次读取128位更改为256位。

int accumulate_i32(__m256i x)
{
    auto tmp1 = _mm256_srli_si256(x,8);
    x = _mm256_add_epi32(x,tmp1);
    auto tmp2 = _mm256_extractf128_si256(x,1);
    tmp2 = _mm_add_epi32(tmp2,_mm256_castsi256_si128(x));
    
    return _mm_cvtsi128_si32(tmp2) + _mm_extract_epi32(tmp2,1);
}

int inner_product_avx(const unsigned char* vec1,unsigned int length)
{    
    constexpr unsigned int memoryAlignmentBytes = 32;
    constexpr unsigned int bytesPerPack = 256 / 8;

    assert((reinterpret_cast<std::uintptr_t>(vec1) % memoryAlignmentBytes) == 0);
    assert((reinterpret_cast<std::uintptr_t>(vec2) % memoryAlignmentBytes) == 0);    
    

    // compute middle part via AVX2    
    unsigned int packCount = length / bytesPerPack;
    const __m256i zeros = _mm256_setzero_si256();
    auto sumlo = _mm256_setzero_si256();
    auto sumhi = _mm256_setzero_si256();

    for (unsigned int packIdx = 0; packIdx < packCount; ++packIdx) {
        auto x1 = _mm256_load_si256((const __m256i*)vec1);
        auto x2 = _mm256_load_si256((const __m256i*)vec2);

        auto x1lo = _mm256_unpacklo_epi8(x1,zeros);
        auto x1hi = _mm256_unpackhi_epi8(x1,zeros);
        auto x2lo = _mm256_unpacklo_epi8(x2,zeros);
        auto x2hi = _mm256_unpackhi_epi8(x2,zeros);
        
        auto tmplo = _mm256_madd_epi16(x1lo,x2lo);
        auto tmphi = _mm256_madd_epi16(x1hi,x2hi);

        sumlo = _mm256_add_epi32(sumlo,tmplo);
        sumhi = _mm256_add_epi32(sumhi,tmphi);

        vec1 += bytesPerPack;
        vec2 += bytesPerPack;
    }

    int sum = accumulate_i32(sumlo) + accumulate_i32(sumhi);

    
    // compute remaining part that cannot be represented as a 
    // whole packed integer    
    unsigned int packRestCount = length % bytesPerPack;
    for (size_t i = packRestCount; i > 0; --i,++vec2)
        sum += int(*vec1) * int(*vec2);

    
    return sum;
}

这大约花费了简单C ++实现的20％的时间（见下文）。考虑到AVX代码可同时处理16个16位整数的事实，我期望会有更高的增益。 -这合理吗？还是我错过了什么？

展开AVX代码中的最后一个循环并不会减少计算时间。

int inner_product_simple(const unsigned char* vec1,size_t length)
{
    int sum = 0;

    for (size_t i = length; i > 0; --i,++vec2)
        sum += int(*vec1) * int(*vec2);

    return sum;
}

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

avxc++correlationsimd