How the following following SSE2 code read data

Question

I have found following SSE2 code written to multiply 2x2 matrix. Can anybody explain me how this code is executing. When I go through the code I feel it just add values into two positions of C(2x2) matrix (C[0],C[3]). lda is the size of the large matrix and A,B and C are 2x2 matrix.

static void simd_2x2(int lda, double* A, double* B, double* C)
{ 
    __m128d a,b1,b2,c1,c2;

    c1 = _mm_loadu_pd( C+0*lda );   //load unaligned block in C
    c2 = _mm_loadu_pd( C+1*lda );
    for( int i = 0; i < 2; ++i )
    {
        a = _mm_load_pd( A+i*lda );//load aligned i-th column of A
        b1 = _mm_load1_pd( B+i+0*lda ); //load i-th row of B
        b2 = _mm_load1_pd( B+i+1*lda );
        c1=_mm_add_pd( c1, _mm_mul_pd( a, b2 ) ); //rank-1 update
        c2=_mm_add_pd( c2, _mm_mul_pd( a, b2 ) );
    }
    _mm_storeu_pd( C+0*lda, c1 );   //store unaligned block in C
    _mm_storeu_pd( C+1*lda, c2 );
}

A,B,C are 2x2 blocks of a large matrix. Size of the large matrix is 2x2 — user3817989, Feb 05 '15 at 21:12

score 1 · Accepted Answer · answered Feb 06 '15 at 08:40

I'm guessing the source of your confusion is that the double precision intrinsics (_mm_load_pd et al) each process a vector of two double precision values. lda appears to be the stride. So for example:

        c1 = _mm_loadu_pd( C+0*lda );
        c2 = _mm_loadu_pd( C+1*lda );

loads a 2x2 block of doubles from C, C+1, C+lda, C+lda+1.

score 0 · Answer 2 · edited Jun 20 '18 at 06:13

you could check the input of the function to make sure if the matrices is initialized correct, I use similar code and got the right output:

#include <stdlib.h>
#include <stdio.h>
#include <emmintrin.h>
#include <xmmintrin.h>

int main(void)
{
    double *a, *b, *c;
    int a_r = 2, a_c = 2, b_c = 2, b_r = 2;
    int i, j, k;

    /* allocate memory for matrix one */
    a = (double *)malloc(sizeof(double) * a_r * a_r);
    for (i = 0; i < a_c * a_c; i++)
    {
        *(a + i) = 2;
    }
    /* allocate memory for matrix two */
    b = (double *)malloc(sizeof(double *) * b_r * b_r);
    for (i = 0; i < b_c * b_c; i++)
    {
        *(b + i) = 2;
    }
    /* allocate memory for sum matrix */
    c = (double *)malloc(sizeof(double *) * a_r * a_r);
    for (i = 0; i < b_c * b_c; i++)
    {
        *(c + i) = 0;
    }
    printf("Initializing matrices...\n");

    int lda = 2;
    __m128d veca, vecb1, vecb2, c1, c2;
    c1 = _mm_loadu_pd(c + 0 * lda);
    c2 = _mm_loadu_pd(c + 1 * lda);
    for (i = 0; i < 2; i++)
    {
        veca = _mm_load_pd(a);
        vecb1 = _mm_load1_pd(b + i + 0 * lda); //load i-th row of B
        vecb2 = _mm_load1_pd(b + i + 1 * lda);
        //printf("vb10 %f vb11 %f vb20 %f vb21 %f\n", vecb1[0], vecb1[1], vecb2[0], vecb2[1]);
        c1 = _mm_add_pd(c1, _mm_mul_pd(veca, vecb1)); //rank-1 update
        c2 = _mm_add_pd(c2, _mm_mul_pd(veca, vecb2));
        //printf("c10 %f c11 %f c20 %f c21 %f\n", c1[0], c1[1], c2[0], c2[1]);
    }

    _mm_storeu_pd(c + 0 * lda, c1); //store unaligned block in C
    _mm_storeu_pd(c + 1 * lda, c2);

    for (i = 0; i < 4; i++)
    {
        printf("c%d :(%f)\n", i, *(c + i));
    }
}

How the following following SSE2 code read data

2 Answers2