使用SSE内在函数的矩阵乘法

Question

我正在尝试使用SSE来进行矩阵乘法。我已经为4x4矩阵编写了一个简单的程序。一切似乎都很好，但是当我打印结果时，它的一些垃圾值。请帮忙弄清楚问题。其次，程序在我释放内存时停止工作，而不是程序的正确结束。

#include <stdlib.h> 
#include <stdio.h> 
#include <time.h> 
#include <float.h> 
#include <xmmintrin.h> 

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

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

    //initializing first matrix 
    for (i = 0; i < a_r; i++) { 
     for (j = 0; j < a_c; j++) { 
      a[i][j] = 2; 
     } 
    } 
    // initializing second matrix 
    for (i = 0; i < b_r; i++) { 
     for (j = 0; j < b_c; j++) { 
      b[i][j] = 2; 
     } 
    } 
    /* initialize product matrix */ 
    for (i = 0; i < a_r; i++) { 
     for (j = 0; j < b_c; j++) { 
      c[i][j] = 0; 
     } 
    } 

    int count = 0; 
    /* multiply matrix one and matrix two */ 
    for (i = 0; i < a_r; i++) { 
     for (j = 0; j < a_c; j++) { 
      count = 0; 
      __m128 result = _mm_setzero_ps(); 
      for (k = 0; k < 4; k += 4) { 
       __m128 row1 = _mm_loadu_ps(&a[i][k]); 
       __m128 row2 = _mm_loadu_ps(&b[k][j]); 
       result = _mm_mul_ps(row1, row2); 

       for (int t = 1; t < 4; t++) { 
        __m128 row3=_mm_loadu_ps(&a[t * 4]); 
        __m128 row4=_mm_loadu_ps(&b[i][t]); 
        __m128 row5 = _mm_mul_ps(row3,row4); 
        result = _mm_add_ps(row5, result); 
       } 
       _mm_storeu_ps(&c[i][j], result); 
      } 
     } 
    } 
    printf("******************************************************\n"); 
    printf ("Done.\n"); 

    for (i = 0; i < a_r ; i++) { 
     for (j = 0; j < b_c; j++) { 
      printf ("%f ", c[i][j]); // issue here when I print results. 
     } 
     printf("\n"); 
    }  // Here program stops working. 

    /*free memory*/ 
    for (i = 0; i < a_r; i++) { 
     free(a[i]); 
    } 
    free(a); 
    for (i = 0; i < a_c; i++) { 
     free(b[i]); 
    } 
    free(b); 
    for (i = 0; i < b_c; i++) { 
     free(c[i]); 
    } 
    free(c); 
} 

请看看输出矩阵打印的地址。如何获得对齐的地址，我有_aligned_malloc，但仍然没有对齐。

Answer 1

用于基质间接指针的分配不正确。应改为：

a = (float **)malloc(sizeof(float*) * a_r); 

写这些分配一个更安全的方法是这样的：

a = malloc(sizeof(*a) * a_r); 

需要注意的是，你可以分配2D直接矩阵：

float (*a)[4][4] = malloc(sizeof(*a)); 

或者更好的，如科迪灰色建议：

float (*a)[4][4] = _aligned_malloc(sizeof(*a)); 

_aligned_malloc是确保SSE操作数正确对齐的非标准函数。

如果事实上你可能甚至不需要与malloc()分配这些矩阵：

float a[4][4]; 

但随着后者的选择，你必须确保对SSE操作成功的正确对齐。

的代码的其余部分有其他问题：

• void main()不正确。它应该是int main(void)

• 第二个矩阵操作数应该转置，以便您可以一次读取多个值。第二次加载将变为：

  __m128 row2 = _mm_loadu_ps(&b[j][k]); 
  

• 求和阶段似乎也不正确。而最终的专卖店肯定是不正确，应该只是：

  c[i][j] = sum;