c - 使用 SSE 对任意大小的输入矩阵和 vector 进行 vector 矩阵和矩阵矩阵乘法

Question

我正在尝试使用 SSE Intrinsic 进行 vector 矩阵乘法以及矩阵矩阵乘法，但如果我尝试执行除 4 的倍数之外的任何操作，我会收到一条错误消息“段错误”。无法弄清楚为什么，它对其他任何东西都不起作用。？请提出更改建议，以便它适用于任何大小的输入。？

以下是我的实现:

#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <xmmintrin.h>
#include <time.h>
#include <omp.h>  

/*****************************************************
the following function generates a "size"-element vector
and a "size x size" matrix
 ****************************************************/
void matrix_vector_gen(int size, float *matrix, float *vector){
    int i;
    for (i = 0; i < size*size; i++){
        vector[i] = i*1.2f + 1;//((float)rand())/65535.0f;
        printf("%f \n ", vector[i]);
    }
    for (i = 0; i < size*size; i++){
        matrix[i] = i*1.3f + 1;//((float)rand())/5307.0f;
        printf("%f \n ", matrix[i]);
    }
}

/****************************************************
the following function calculate the below equation
   vector_out = vector_in x matrix_in
 ***************************************************/
void matrix_mult_sq(int size, float *vector_in,
               float *matrix_in, float *vector_out){
    int i, j, k;
    for (i = 0; i < size; i++)
    {
        for (j = 0; j < size; j++)
        {
            vector_out[size*i + j] = 0.0;
            for (k = 0; k < size; k++)
                vector_out[size*i + j] += vector_in[size*i + k] * matrix_in[size*k + j];
        }
    }
}

void matrix_mult_sse(int size, float *vector_in,
    float *matrix_in, float *vector_out){
    __m128 a_line, b_line, r_line;
    int i, j, k, l;
    for (k = 0; k < size; k++)
    {

        for (i = 0; i < size; i += 4){
            j = 0;
            b_line = _mm_load_ps(&matrix_in[i]); // b_line = vec4(matrix[i][0])
            a_line = _mm_set1_ps(vector_in[j + k*size]);      // a_line = vec4(vector_in[0])
            r_line = _mm_mul_ps(a_line, b_line); // r_line = a_line * b_line
            for (j = 1; j < size; j++) {
                b_line = _mm_load_ps(&matrix_in[j*size + i]); // a_line = vec4(column(a, j))
                a_line = _mm_set1_ps(vector_in[j + k*size]);  // b_line = vec4(b[i][j])
                // r_line += a_line * b_line
                r_line = _mm_add_ps(_mm_mul_ps(a_line, b_line), r_line);
            }
            _mm_store_ps(&vector_out[i + k*size], r_line);     // r[i] = r_line
        }
    }
    for (l=0; l < size*size; l++)
    {
        printf("%f \n", vector_out[l]);
    }
}

int main(int argc, char *argv[]){
  if(argc < 2){
    printf("Usage: %s matrix/vector_size\n", argv[0]);
    return 0;
  }

  int size = atoi(argv[1]);
  if(size%4 != 0){
    printf("This version implements for ""size = 4*n"" only\n");
    return 0;
  }

  float *vector = (float *)memalign(sizeof(float)*4, sizeof(float)*size);//(float *)malloc(sizeof(float)*size);
  if(vector==NULL){
    printf("can't allocate the required memory for vector\n");
    return 0;
  }

  float *matrix = (float *)memalign(sizeof(float)*4, sizeof(float)*size*size);
  if(matrix==NULL){
    printf("can't allocate the required memory for matrix\n");
    free(vector);
    return 0;
  }

  float *result_sq = (float *)memalign(sizeof(float)*4, sizeof(float)*size);
  if(result_sq==NULL){
    printf("can't allocate the required memory for result_sq\n");
    free(vector);
    free(matrix);
    return 0;
  }

  float *result_pl = (float *)memalign(sizeof(float)*4, sizeof(float)*size);
  if(result_pl==NULL){
    printf("can't allocate the required memory for result_pl\n");
    free(vector);
    free(matrix);
    free(result_sq);
    return 0;
  }

  matrix_vector_gen(size, matrix, vector);

  double time_sq;
  double time_sse;

  time_sq = omp_get_wtime();
  matrix_mult_sq(size, vector, matrix, result_sq);
  time_sq = omp_get_wtime() - time_sq;

  time_sse = omp_get_wtime();
  matrix_mult_sse(size, vector, matrix, result_pl);
  time_sse = omp_get_wtime() - time_sse;

  printf("SEQUENTIAL EXECUTION: %f (sec)\n",time_sq);
  printf("PARALLEL EXECUTION: %f (sec)\n", time_sse);

  //check
  /*int i;
  for(i=0; i<size; i++)
    if((int)result_sq[i] != (int)result_pl[i]){
      printf("wrong at position %d\n", i);
      free(vector);
      free(matrix);
      free(result_sq);
      free(result_pl);
      return 0;
    }*/

  free(vector);
  free(matrix);
  free(result_sq);
  free(result_pl);
  return 1;
}

Answer 1

看来您专门使用 mm_load_ps 和 mm_store_ps 来加载和存储，它们在一条指令中加载和存储 4 个 float 。

由于您的容器(矩阵和 vector )的大小不一定是 4 个浮点(16 字节)的倍数，因此这是不正确的。

memalign 确保指针对齐(此处为 16 字节)，但不会在末尾保留填充，以便分配的 block 大小是 16 字节的倍数。

例如，存储一个 5 维 vector 时，该 vector 在内存中只分配了 20 个字节，但你写入了 32 个字节(两次 mm_store_ps 操作)

此外，这似乎是不正确的:

_mm_store_ps(&vector_out[i + k*size], r_line);

如果我是对的，你想在这里存储一个 float 。不是四个打包的 float 。