c++ - OpenMP C++ 矩阵乘法

Question

我想并行化以下代码。尤其是这些 for 循环，因为它是最昂贵的操作。

      for (i = 0; i < d1; i++)
         for (j = 0; j < d3; j++)
             for (k = 0; k < d2; k++)
             C[i][j] = C[i][j] + A[i][k] * B[k][j];  

这是我第一次尝试使用 OpenMP 并行化代码。我已经尝试了几种方法，但我总是以比使用串行版本更糟糕的运行时间告终。 如果您能告诉我代码或编译指示是否有问题，那就太好了...

      #include <omp.h>
      #include <stdio.h>
      #include <stdlib.h>
      //#include <stdint.h>

      // ---------------------------------------------------------------------------
      // allocate space for empty matrix A[row][col]
      // access to matrix elements possible with:
      // - A[row][col]
      // - A[0][row*col]


      float **alloc_mat(int row, int col)
      {
          float **A1, *A2;

          A1 = (float **)calloc(row, sizeof(float *));      // pointer on rows
          A2 = (float *)calloc(row*col, sizeof(float));    // all matrix elements

          //#pragma omp parallel for
          for (int i=0; i<row; i++)
              A1[i] = A2 + i*col;

          return A1;
      }

      // ---------------------------------------------------------------------------
      // random initialisation of matrix with values [0..9]

      void init_mat(float **A, int row, int col)
      {   
          //#pragma omp parallel for
          for (int i = 0; i < row*col; i++)
              A[0][i] = (float)(rand() % 10);
      }

      // ---------------------------------------------------------------------------
      // DEBUG FUNCTION: printout of all matrix elements

      void print_mat(float **A, int row, int col, char *tag)
      {
          int i, j;

          printf("Matrix %s:\n", tag);
          for (i = 0; i < row; i++)
          {
              //#pragma omp parallel for
              for (j=0; j<col; j++) 
                  printf("%6.1f   ", A[i][j]);
              printf("\n"); 
          }
      }

      // ---------------------------------------------------------------------------

      int main(int argc, char *argv[])
      {
          float **A, **B, **C;  // matrices
          int d1, d2, d3;         // dimensions of matrices
          int i, j, k;          // loop variables


          double start, end;
          start = omp_get_wtime();

          /* print user instruction */
          if (argc != 4)
          {
              printf ("Matrix multiplication: C = A x B\n");
              printf ("Usage: %s <NumRowA>; <NumColA> <NumColB>\n",argv[0]); 
               return 0;
           }

           /* read user input */
           d1 = atoi(argv[1]);      // rows of A and C
           d2 = atoi(argv[2]);     // cols of A and rows of B
           d3 = atoi(argv[3]);     // cols of B and C

           printf("Matrix sizes C[%d][%d] = A[%d][%d] x B[%d][%d]\n", 
           d1, d3, d1, d2, d2, d3);

           /* prepare matrices */
           A = alloc_mat(d1, d2);
           init_mat(A, d1, d2); 
           B = alloc_mat(d2, d3);
           init_mat(B, d2, d3);
           C = alloc_mat(d1, d3);   // no initialisation of C, 
       //because it gets filled by matmult

           /* serial version of matmult */
           printf("Perform matrix multiplication...\n");



           int sum;
           //#pragma omp parallel
           //{
               #pragma omp parallel for collapse(3) schedule(guided)
               for (i = 0; i < d1; i++)
                   for (j = 0; j < d3; j++)
                       for (k = 0; k < d2; k++){
                       C[i][j] = C[i][j] + A[i][k] * B[k][j];
                       }
           //}


           end = omp_get_wtime();


           /* test output */
           print_mat(A, d1, d2, "A"); 
           print_mat(B, d2, d3, "B"); 
           print_mat(C, d1, d3, "C"); 

           printf("This task took %f seconds\n", end-start);
           printf ("\nDone.\n");

           return 0;
       }

Answer 1

正如@genisage 在评论中建议的那样，矩阵的大小可能足够小，以至于初始化额外线程的开销大于通过并行计算矩阵乘法所节省的时间。但是，请考虑下图，其中包含我通过使用和不使用 OpenMP 运行您的代码获得的数据。

我使用了从 n=10 到 n=1000 的方阵。请注意，在 n=50 和 n=100 之间的某处，并行版本如何变得更快。

然而，在尝试编写快速矩阵乘法时，还有其他问题需要考虑，这主要与有效使用缓存有关。首先，您连续分配整个矩阵(这很好)，但仍然使用两个指针重定向来访问数据，这是不必要的。此外，您的矩阵以行主要格式存储，这意味着您正在连续访问 A 和 C 中的数据，而不是 B 中的数据。不是显式存储 B 并将 A 的行与 B 的列相乘，您会得到通过存储 B 转置并将 A 的一行元素与 B 转置的一行相乘来更快地乘法。

这是一个仅针对 A*B 的优化，然而，您的代码中可能还有其他地方存储 B 比 B 转置更好，在这种情况下通常执行 matrix multiplication by blocking可以导致更好的缓存利用率