下面的SSE2代码如何读取数据

Question

我发现下面的SSE2代码写成乘以2x2矩阵。任何人都可以解释我的代码是如何执行的。当我通过代码时，我感觉它只是将值添加到C（2x2）矩阵（C [0]，C [3]）的两个位置中。“lda是大矩阵的大小，A，B和C是2×2矩阵”

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); 
} 

Answer 1

我猜你的混乱的根源在于，双精度内在（_mm_load_pd等），每个处理两个双精度值的向量。似乎是迈步。例如：

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

从C，C + 1，C + Ida，C + Ida + 1加载一个2x2块双打。

Answer 2

你可以检查函数的输入，以确保如果矩阵被初始化正确的，我用类似的代码，并得到了正确的输出：

#include <stdlib.h> 
#include <stdio.h> 
#include <time.h> 
#include <float.h> 
#include <emmintrin.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)); 
} 
}