Optimizing code using Intel SSE intrinsics for vectorization

Question

This is my very first time working with SSE intrinsics. I am trying to convert a simple piece of code into a faster version using Intel SSE intrinsic (up to SSE4.2). I seem to encounter a number of errors.

The scalar version of the code is: (simple matrix multiplication)

     void mm(int n, double *A, double *B, double *C)
     {
        int i,j,k;
        double tmp;

        for(i = 0; i < n; i++)
            for(j = 0; j < n; j++) {
                    tmp = 0.0;
                    for(k = 0; k < n; k++)
                            tmp += A[n*i+k] *
                                   B[n*k+j];
                    C[n*i+j] = tmp;

              }
            }

This is my version: I have included #include <ia32intrin.h>

      void mm_sse(int n, double *A, double *B, double *C)
      {
        int i,j,k;
        double tmp;
        __m128d a_i, b_i, c_i;

        for(i = 0; i < n; i++)
            for(j = 0; j < n; j++) {
                    tmp = 0.0;
                    for(k = 0; k < n; k+=4)
                            a_i = __mm_load_ps(&A[n*i+k]);
                            b_i = __mm_load_ps(&B[n*k+j]);
                            c_i = __mm_load_ps(&C[n*i+j]);

                            __m128d tmp1 = __mm_mul_ps(a_i,b_i);
                            __m128d tmp2 = __mm_hadd_ps(tmp1,tmp1);
                            __m128d tmp3 = __mm_add_ps(tmp2,tmp3);
                            __mm_store_ps(&C[n*i+j], tmp3);

            }
         }

Where am I going wrong with this? I am getting several errors like this:

mm_vec.c(84): error: a value of type "int" cannot be assigned to an entity of type "__m128d" a_i = __mm_load_ps(&A[n*i+k]);

This is how I am compiling: icc -O2 mm_vec.c -o vec

Can someone please assist me converting this code accurately. Thanks!

UPDATE:

According to your suggestions, I have made the following changes:

       void mm_sse(int n, float *A, float *B, float *C)
       {
         int i,j,k;
         float tmp;
         __m128 a_i, b_i, c_i;

         for(i = 0; i < n; i++)
            for(j = 0; j < n; j++) {
                    tmp = 0.0;
                    for(k = 0; k < n; k+=4)
                            a_i = _mm_load_ps(&A[n*i+k]);
                            b_i = _mm_load_ps(&B[n*k+j]);
                            c_i = _mm_load_ps(&C[n*i+j]);

                            __m128 tmp1 = _mm_mul_ps(a_i,b_i);
                            __m128 tmp2 = _mm_hadd_ps(tmp1,tmp1);
                            __m128 tmp3 = _mm_add_ps(tmp2,tmp3);
                            _mm_store_ps(&C[n*i+j], tmp3);


            }
        }

But now I seem to be getting a Segmentation fault. I know this perhaps because I am not accessing the array subscripts properly for array A,B,C. I am very new to this and not sure how to proceed with this.

Please help me determine the correct approach towards handling this code.

Eleanoraeleanore

__m128d tmp3 = __mm_add_ps(tmp2,tmp3); is total nonsense: you're using the variable as it's declared, thus uninitialized. See Get sum of values stored in __m256d with SSE/AVX for correct + efficient horizontal sums. (But really, you don't want to do that in the inner loop. Apply SIMD over the 2nd loop, for example, to calculate for row*column dot products in parallel. Also, _ps is packed-single; not what you want for double*. – Alegre

Answer 1

The error you're seeing is because you have too many underscores in the function names, e.g.:

__mm_mul_ps

should be:

_mm_mul_ps // Just one underscore up front

so the C compiler is assuming they return int since it hasn't seen a declaration.

Beyond this though there's further problems - you seem to be mixing calls to double and single float variants of the same instruction.

For example you have:

__m128d a_i, b_i, c_i;

but you call:

__mm_load_ps(&A[n*i+k]);

which returns a __m128 not a __m128d - you wanted to call:

_mm_load_pd

instead. Likewise for the other instructions if you want them to work on pairs of doubles.

---

If you're seeing unexplained segmentation faults and in SSE code I'd be inclined to guess that you've got memory alignment problems - pointers passed to SSE intrinsics (mostly¹) need to be 16 byte aligned. You can check this with a simple assert in your code, or check it in a debugger (you expect the last digit of the pointer to be 0 if it's aligned properly).

If it isn't aligned right you need to make sure it is. For things not allocated with new/malloc() you can do this with a compiler extension (e.g. with gcc):

float a[16] __attribute__ ((aligned (16)));

Provided your version of gcc has a max alignment large enough to support this and a few other caveats about stack alignment. For dynamically allocated storage you'll want to use a platform specific extension, e.g. posix_memalign to allocate suitable storage:

float *a=NULL;
posix_memalign(&a, __alignof__(__m128), sizeof(float)*16);

(I think there might be nicer, portable ways of doing this with C++11 but I'm not 100% sure on that yet).

¹ There are some instructions which allow you do to unaligned loads and stores, but they're terribly slow compared to aligned loads and worth avoiding if at all possible.

Dimissory

I am working with icc not gcc. Do you think that handling it like this: a_i = _mm_load_ps(&A[n*i+k]);is the correct approach? The examples I see posted elsewhere (an even on the Intel Intrinsic documentation) have very basic examples. arrays A,B,C have all been allocated with malloc. – Eleanoraeleanore

@Hello_PG The load itself isn't directly wrong. You don't need to load c_i though. For the most part ICC has the same extensions as gcc - I think that's the case for the alignment one, I'm more familiar with GCC than ICC personally though so I qualified it with that and linked to the docs I knew how to find. malloc doesn't guarantee suitable alignment on all platforms, so posix_memalign is probably needed. Did the assert I suggested fail? – Dimissory

When I try to allocate memory for A like this: A = (float*)_aligned_malloc(dimensiondimensionsizeof(float),16); I get a compile error:undefined reference to `aligned_malloc' with icc. this is how i am compiling:icc -O2 mm_vec.c -o vec2 – Eleanoraeleanore

@Hello_PG - it's either aligned_mallloc (note: no underscore) although probably you want to use posix_memalign, not aligned malloc anyway. (You probably want to put compiler warnings on and note when it says about implicit declarations of functions). – Dimissory

Answer 2

You need to make sure that your loads and stores are always accessing 16 byte aligned addresses. Alternatively, if you can't guarantee this for some reason, then use _mm_loadu_ps/_mm_storeu_ps instead of _mm_load_ps/_mm_store_ps - this will be less efficient but will not crash on misaligned addresses.

Sequel