15

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

Peter Cordes
  • 328,167
  • 45
  • 605
  • 847
gpuguy
  • 4,607
  • 17
  • 67
  • 125

3 Answers3

11

You can abuse PSADBW to calculate horizontal sums of bytes without overflow. For example:

pxor    xmm0, xmm0
psadbw  xmm0, [a + 0]     ; sum in 2x 64-bit chunks
pxor    xmm1, xmm1
psadbw  xmm1, [a + 16]
paddw   xmm0, xmm1        ; accumulate vertically
pshufd  xmm1, xmm0, 2     ; bring down the high half
paddw   xmm0, xmm1   ; low word in xmm0 is the total sum
; movd  eax, xmm0    ; higher bytes are zero so efficient dword extract is fine

Intrinsics version:

#include <immintrin.h>
#include <stdint.h>

// use loadu instead of load if 16-byte alignment of a[] isn't guaranteed
unsigned sum_32x8(const uint8_t a[32])
{
    __m128i zero = _mm_setzero_si128();
    __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_epi32(sum0, sum1);
    __m128i totalsum = _mm_add_epi32(sum2, _mm_shuffle_epi32(sum2, 2));
    return _mm_cvtsi128_si32(totalsum);
}

This portably compiles back to the same asm, as you can see on Godbolt.

The reinterpret_cast<const __m128i*> is necessary because Intel intrinsics before AVX-512 for integer vector load/store take __m128i* pointer args, instead of a more convenient void*. Some prefer more compact C-style casts like _mm_loadu_si128( (const __m128*) &a[16] ) as a style choice.

16 vs. 32 vs. 64-bit SIMD element size doesn't matter much; 16 and 32 are equally efficient on all machines, and 32-bit will avoid overflow even if you use this for summing much larger arrays. (paddq is slower on some old CPUs like Core 2; https://agner.org/optimize/ and https://uops.info/) Extracting as 32-bit is definitely more efficient than _mm_extract_epi16 (pextrw).

Peter Cordes
  • 328,167
  • 45
  • 605
  • 847
harold
  • 61,398
  • 6
  • 86
  • 164
  • Could you please write the Intel® C++ Compiler Intrinsic Equivalents for the above? – gpuguy Jun 12 '12 at 03:07
  • To use this same trick for `int8_t` (instead of `uint8_t`): range shift to unsigned (xor with 0x80), then subtract `16 * 0x80` from the total. See [this patch for Agner Fog's vector class library for an example with intrinsics](https://github.com/pcordes/vectorclass/commit/630ca802bb1abefd096907f8457d090c28c8327b). The same idea works for [an AVX2 ymm vector](https://github.com/pcordes/vectorclass/commit/11aa77071d25d7d93090789006250f8992f44272)). – Peter Cordes Jun 11 '16 at 08:36
5

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.

Paul R
  • 208,748
  • 37
  • 389
  • 560
  • How do I make sure that a[] is 16 byte aligned? In SSE is there any thing similar to __align__(16) in CUDA? – gpuguy Jun 08 '12 at 03:01
  • It depends what compiler and OS you are using - e.g. for gcc with non-dynamic allocations use `__attribute__ ((aligned(16)))` - for dynamic allocations on Linux use `memalign()` or `posix_memalign()`. – Paul R Jun 08 '12 at 06:26
  • going to have to downvote this; it works but `psadbw` is the Right Answer. For signed `int8_t`, you can range-shift to unsigned with `xor` to add 0x80 to each byte, and subtract out `16 * 0x80` from the result. (See [this patch for Agner Fog's vector class library](https://github.com/pcordes/vectorclass/commit/630ca802bb1abefd096907f8457d090c28c8327b) for example with intrinsics. The same idea works for [an AVX2 ymm vector](https://github.com/pcordes/vectorclass/commit/11aa77071d25d7d93090789006250f8992f44272)). But the OP here seems to already have unsigned, so psadbw is a huge win. – Peter Cordes Jun 11 '16 at 00:26
  • 1
    @PeterCordes: well I did say it was "a bit long-winded". ;-) But yes, harold's answer is the better solution (and I've upvoted it of course). I should probably delete this one as it doesn't really serve any useful purpose. – Paul R Jun 11 '16 at 07:14
0

There is one more way to find the sum of all elements of an array using SSE instructions. The code uses the following SSE constructs.

  • __m256 register
  • _mm256_store_ps(float *a, __m256 b)
  • _mm256_add_ps(__m256 a, __m256 b)

The code works for any sized array of floats.

float sse_array_sum(float *a, int size)
{
    /*
     *   sum += a[i] (for all i in domain)
     */

    float *sse_sum, sum=0;
    if(size >= 8)
    {
        // sse_sum[8]
        posix_memalign((void **)&sse_sum, 32, 8*sizeof(float));

        __m256 temp_sum;
        __m256* ptr_a = (__m256*)a;
        int itrs = size/8-1;

        // sse_sum[0:7] = a[0:7]
        temp_sum = *ptr_a;
        a += 8;
        ptr_a++;

        for(int i=0; i<itrs; i++, ptr_a++, a+=8)
            temp_sum = _mm256_add_ps(temp_sum, *ptr_a);

        _mm256_store_ps(sse_sum, temp_sum);
        for(int i=0; i<8; i++)  sum += sse_sum[i];
    }

    // if size is not divisible by 8
    int rmd_itrs = size%8;
    // Note: a is pointing to remainder elements
    for(int i=0; i<rmd_itrs; i++)   sum += a[i];

    return sum;
}


float seq_array_sum(float *a, int size)
{
    /*
     *  sum += a[i] (for all i)
     */

    float sum = 0;
    for(int i=0; i<size; i++)   sum += a[i];
    return sum;
}

Benchmark:

size = 64000000
a[i] = 3141592.65358 for all i in domain

sequential version time: 194ms
SSE version time: 49ms

Machine specification:

Thread(s) per core: 2
Core(s) per socket: 2
Socket(s): 1
CPU MHz: 1700.072
OS: Ubuntu

dwijesh522
  • 29
  • 4
  • 1
    First of all, this question is about summing `uint8_t[]`. More importantly, dynamically allocating temporary storage is definitely not a win. Just use a `__m256 sum` temporary like a normal person. (Or better, use multiple accumulators to hide FP latency. If you want to use a local "array" of them, compilers will typically optimize that away and keep them all in registers). And definitely don't `_mm256_store_ps()` inside the inner loop, although that might optimize away, too. – Peter Cordes Mar 31 '20 at 06:14
  • 1
    Also, your code isn't safe if `a` isn't aligned by 32; you dereference a `__m256*` instead of using `_mm256_loadu_ps`. – Peter Cordes Mar 31 '20 at 06:14
  • 1
    An AVX `__m256` version of [SSE reduction of float vector](https://stackoverflow.com/q/17761154) is what you want. Note how simple its inner loop is. For the horizontal sum at the end, see [Fastest way to do horizontal SSE vector sum (or other reduction)](https://stackoverflow.com/q/6996764). Storing to a `alignas(32) float tmp[8]` would be one option, but you can do better with shuffles. Also, unroll the main loop with multiple accumulators to hide FP latency like in [this Q&A](https://stackoverflow.com/q/45113527) – Peter Cordes Mar 31 '20 at 06:25
  • Thanks for the feed-back. You are right. Next time I will take care of 1) being specific to answer 2) getting the proper background. Actually I just started using SSE instructions and thought this might help someone. – dwijesh522 Mar 31 '20 at 09:18
  • 1
    Just FYI: `__m256` stuff requires AVX. It's not totally wrong to call all x86 SIMD stuff in general "SSE", but usually people at least separate SSE* vs. AVX1 / AVX2 / FMA vs. AVX512. Anyway, thanks for trying to help out, but unfortunately this answer doesn't demonstrate good ways to do things, and I'd have to downvote it even if you posted it on [SSE reduction of float vector](https://stackoverflow.com/q/17761154) where it would at least be answering the right question. – Peter Cordes Mar 31 '20 at 09:41
  • 1
    That's normal when you're first learning something, don't feel bad that your first attempt was overcomplicated and unsafe (only works for aligned `a`, and leaks dynamically allocated memory - like I said just use a local array or use shuffles to hsum). But don't that doesn't change the fact that it's not a good example for other beginners or anyone and thus not a good Stack Overflow answer. Perhaps https://codereview.stackexchange.com/ would be a good place to post attempts like this to get feedback on how to do it properly. – Peter Cordes Mar 31 '20 at 09:43