Sum reduction of unsigned bytes without overflow,

Question

I am trying to find sum reduction of 32 elements (each 1 byte data) on an Intel i3 processor. I did this:

s=0; 
for (i=0; i<32; i++)
{
    s = s + a[i];
}  

However, its taking more time, since my application is a real-time application requiring much lesser time. Please note that the final sum could be more than 255.

Is there a way I can implement this using low level SIMD SSE2 instructions? Unfortunately I have never used SSE. I tried searching for sse2 function for this purpose, but it is also not available. Is it (sse) guaranteed to reduce the computation time for such a small-sized problems?

Any suggestions??

Note: I have implemented the similar algorithms using OpenCL and CUDA and that worked great but only when the problem size was big. For small sized problems the cost of overhead was more. Not sure how it works on SSE

Answer 1

This is a bit long-winded but it should still be at least 2x faster than the scalar code:

uint16_t sum_32(const uint8_t a[32])
{
    const __m128i vk0 = _mm_set1_epi8(0);   // constant vector of all 0s for use with _mm_unpacklo_epi8/_mm_unpackhi_epi8
    __m128i v = _mm_load_si128(a);          // load first vector of 8 bit values
    __m128i vl = _mm_unpacklo_epi8(v, vk0); // unpack to two vectors of 16 bit values
    __m128i vh = _mm_unpackhi_epi8(v, vk0);
    __m128i vsum = _mm_add_epi16(vl, vh);
    v = _mm_load_si128(&a[16]);             // load second vector of 8 bit values
    vl = _mm_unpacklo_epi8(v, vk0);         // unpack to two vectors of 16 bit values
    vh = _mm_unpackhi_epi8(v, vk0);
    vsum = _mm_add_epi16(vsum, vl);
    vsum = _mm_add_epi16(vsum, vh);
    // horizontal sum
    vsum = _mm_add_epi16(vsum, _mm_srli_si128(vsum, 8));
    vsum = _mm_add_epi16(vsum, _mm_srli_si128(vsum, 4));
    vsum = _mm_add_epi16(vsum, _mm_srli_si128(vsum, 2));
    return _mm_extract_epi16(vsum, 0);
}

Note that a[] needs to be 16 byte aligned.

You can probably improve on the above code using _mm_hadd_epi16.

Answer 2

You can abuse PSADBW to calculate small horizontal sums quickly.

Something like this: (not tested)

pxor xmm0, xmm0
psadbw xmm0, [a + 0]
pxor xmm1, xmm1
psadbw xmm1, [a + 16]
paddw xmm0, xmm1
pshufd xmm1, xmm0, 2
paddw xmm0, xmm1 ; low word in xmm0 is the total sum

---

Attempted intrinsics version:

I never use intrinsics so this code probably makes no sense whatsoever. The disassembly looked OK though.

uint16_t sum_32(const uint8_t a[32])
{
    __m128i zero = _mm_xor_si128(zero, zero);
    __m128i sum0 = _mm_sad_epu8(
                        zero,
                        _mm_load_si128(reinterpret_cast<const __m128i*>(a)));
    __m128i sum1 = _mm_sad_epu8(
                        zero,
                        _mm_load_si128(reinterpret_cast<const __m128i*>(&a[16])));
    __m128i sum2 = _mm_add_epi16(sum0, sum1);
    __m128i totalsum = _mm_add_epi16(sum2, _mm_shuffle_epi32(sum2, 2));
    return totalsum.m128i_u16[0];
}