sseを使用したドットプロダクション

Question

#define Size 50000 

void main() 
{ 

    unsigned char *arry1 = (unsigned char*)malloc(sizeof(unsigned char)* Size); 
    unsigned char *arry2 = (unsigned char*)malloc(sizeof(unsigned char)* Size); 
    unsigned int *result = (unsigned int*)malloc(sizeof(unsigned int)* Size); 


    for (int i = 0; i < 16; i++) 
    { 
     arry1[i] = i; 
     arry2[i] = i; 
    } 


    __m128i Z = _mm_setzero_si128(); 
    __m128i vsum = _mm_setzero_si128(); 
    //__m128i dummy = _mm_setzero_si128(); 

    for (int j = 0; j < 16; j += 16) 
    { 
     //printf("%d\n\n", j); 

     __m128i test1 = _mm_setzero_si128(); 
      test1 = _mm_loadu_si128((__m128i*)&arry1[j]); 
      __m128i test2 = _mm_setzero_si128(); 
      test2 = _mm_loadu_si128((__m128i*)&arry2[j]); 

     __m128i s16L = _mm_unpacklo_epi8(test1, Z); 
     __m128i s16H = _mm_unpackhi_epi8(test1, Z); 
     __m128i s32LL = _mm_unpacklo_epi16(s16L, Z); 
     __m128i s32LH = _mm_unpackhi_epi16(s16L, Z); 
     __m128i s32HL = _mm_unpacklo_epi16(s16H, Z); 
     __m128i s32HH = _mm_unpackhi_epi16(s16H, Z); 

     __m128i t16L = _mm_unpacklo_epi8(test2, Z); 
     __m128i t16H = _mm_unpackhi_epi8(test2, Z); 
     __m128i t32LL = _mm_unpacklo_epi16(t16L, Z); 
     __m128i t32LH = _mm_unpackhi_epi16(t16L, Z); 
     __m128i t32HL = _mm_unpacklo_epi16(t16H, Z); 
     __m128i t32HH = _mm_unpackhi_epi16(t16H, Z); 

     __m128 s1 = _mm_cvtepi32_ps(s32LL); 
     __m128 s2 = _mm_cvtepi32_ps(s32LH); 
     __m128 s3 = _mm_cvtepi32_ps(s32HL); 
     __m128 s4 = _mm_cvtepi32_ps(s32HH); 

     __m128 t1 = _mm_cvtepi32_ps(t32LL); 
     __m128 t2 = _mm_cvtepi32_ps(t32LH); 
     __m128 t3 = _mm_cvtepi32_ps(t32HL); 
     __m128 t4 = _mm_cvtepi32_ps(t32HH); 

     s1 = _mm_mul_ps(s1, t1); 
     s2 = _mm_mul_ps(s2, t2); 
     s3 = _mm_mul_ps(s3, t3); 
     s4 = _mm_mul_ps(s4, t4); 

     s1 = _mm_hadd_ps(s1, s2);//41,13 
     s3 = _mm_hadd_ps(s3, s4); //313,221 

     vsum = _mm_cvtps_epi32(s3); 

     for (int k = 0; k < 16; k++) 
     { 
      printf("%u\n", (unsigned char)vsum.m128i_i8[k]); 
     } 

     s1 = _mm_hadd_ps(s1, s3); //734, 14 
     s1 = _mm_hadd_ps(s1, s1); //1100,140 
     s1 = _mm_hadd_ps(s1, s1); //1240 


    } 


} 

私はsseを使ってドットプロダクションを進めています。私は_mm_mul_psと_mm_hadd_ps命令ではなく_mm_dp_psを使用しています。 _mm_hadd_ps機能後の値が255を超えると、間違った値が表示されます。

例えば、s3の正しい値は{0,0,0,421,0,0,0,313,0,0,0,221,0,0,0,145}です。 ただし、{0,0,1,165,0,0,1,57,0,0,0,221,0,0,0,145}が印刷されます。これは私がarry1、arry2と言った結果ですか？unsigned char？私は255が8ビットで最大値であることを知っています。

Answer 1

私はここにいくつかの問題を見ることができます：あなたは50000のuint8_t値の内積を計算したい場合は

1）、それはOKです。しかし、70000値のドット積は、uint32_t型のオーバーフローを引き起こす可能性があります。したがって、uint64_tを使用する方が良い解決策です。

2）整数ベクトルの内積を計算するには、浮動小数点数を使用する必要はありません。計算には整数だけを使用する方が効果的です。

2つのuint8_tのベクトルのドット積を計算SSE2関数の例がある：

#include <algorithm> 
#include <emmintrin.h> 

const __m128i Z = _mm_setzero_si128(); 
const size_t A = sizeof(__m128i); 
const size_t B = 0x10000; 

inline __m128i DotProduct32(const uint8_t * a, const uint8_t * b) 
{ 
    __m128i _a = _mm_loadu_si128((__m128i*)a); 
    __m128i _b = _mm_loadu_si128((__m128i*)b); 
    __m128i lo = _mm_madd_epi16(_mm_unpacklo_epi8(_a, Z), _mm_unpacklo_epi8(_b, Z)); 
    __m128i hi = _mm_madd_epi16(_mm_unpackhi_epi8(_a, Z), _mm_unpackhi_epi8(_b, Z)); 
    return _mm_add_epi32(lo, hi); 
} 

inline __m128i HorizontalSum32(__m128i a) 
{ 
    return _mm_add_epi64(_mm_unpacklo_epi32(a, Z), _mm_unpackhi_epi32(a, Z)); 
} 

inline uint64_t ExtractSum64(__m128i a) 
{ 
    uint64_t _a[2]; 
    _mm_storeu_si128((__m128i*)_a, a); 
    return _a[0] + _a[1]; 
} 

void DotProduct(const uint8_t * a, const uint8_t * b, size_t size, uint64_t * sum) 
{ 
    size_t blockNumber = (size + B - 1)/B; 
    size_t alignedSize = size/A*A; 
    size_t i = 0; 

    __m128i sum64 = Z; 
    for (size_t block = 0; block < blockNumber; ++i) 
    { 
     __m128i sum32 = Z; 
     for (size_t blockEnd = std::min(alignedSize, i + B); i < blockEnd; i += A) 
      sum32 = _mm_add_epi32(sum32, DotProduct32(a + i, b + i)); 
     sum64 = _mm_add_epi64(sum64, HorizontalSum32(sum32)); 
    } 

    *sum = ExtractSum64(sum64); 
    for (; i < size; ++i) 
     *sum += a[i] * b[i]; 
}