c++ - 矩阵乘法、KIJ 顺序、并行版本比非并行版本慢

Question

我有一项关于并行编程的学校任务，我遇到了很多问题。 我的任务是创建给定矩阵乘法代码的并行版本并测试其性能(是的，它必须按 KIJ 顺序):

void multiply_matrices_KIJ()
{
    for (int k = 0; k < SIZE; k++)
        for (int i = 0; i < SIZE; i++)
            for (int j = 0; j < SIZE; j++)
                matrix_r[i][j] += matrix_a[i][k] * matrix_b[k][j];
}

这是我到目前为止想出的:

void multiply_matrices_KIJ()
{
    for (int k = 0; k < SIZE; k++)
#pragma omp parallel
    {
#pragma omp for schedule(static, 16)
        for (int i = 0; i < SIZE; i++)
            for (int j = 0; j < SIZE; j++)
                matrix_r[i][j] += matrix_a[i][k] * matrix_b[k][j];
    }
}

这就是我发现让我感到困惑的地方。这个代码的并行版本比非并行版本慢 50% 左右。根据矩阵大小，速度差异仅略有不同(测试的 SIZE = 128、256、512、1024、2048 和各种计划版本 - 动态、静态、完全没有它等等)。

有人可以帮助我了解我做错了什么吗？是否可能是因为我使用的是 KIJ 命令，使用 openMP 不会变得更快？

编辑:

我在 Windows 7 PC 上工作，使用 Visual Studio 2015 Community Edition，在 Release x86 模式下编译(x64 也无济于事)。我的 CPU 是:Intel Core i5-2520M CPU @ 2,50GHZ(是的，是的，这是一台笔记本电脑，但我在家用 I7 PC 上得到了相同的结果)

我正在使用全局数组:

float matrix_a[SIZE][SIZE];    
float matrix_b[SIZE][SIZE];    
float matrix_r[SIZE][SIZE];

我正在为矩阵 a 和 b 分配随机(浮点)值，矩阵 r 填充为 0。

到目前为止，我已经用各种矩阵大小(128、256、512、1024、2048 等)测试了代码。对于其中一些，它的目的是不适合缓存。 我当前的代码版本如下所示:

void multiply_matrices_KIJ()
{
#pragma omp parallel 
    {
    for (int k = 0; k < SIZE; k++) {
#pragma omp for schedule(dynamic, 16) nowait
        for (int i = 0; i < SIZE; i++) {
            for (int j = 0; j < SIZE; j++) {
                matrix_r[i][j] += matrix_a[i][k] * matrix_b[k][j];
            }
        }
    }
    }
}

需要说明的是，我知道通过不同的循环顺序我可以获得更好的结果，但事实就是如此——我必须使用 KIJ 顺序。我的任务是并行执行 KIJ for 循环并检查性能提升情况。我的问题是我希望(ed)至少执行得更快一点(比我现在得到的最多快 5-10%)即使它是并行的 I 循环(不能用K 循环，因为我会得到不正确的结果，因为它是 matrix_r[i][j])。

这些是我在使用上面显示的代码时得到的结果(我正在计算数百次并获得平均时间):

大小 = 128

• 序列号:0,000608s
• 并行 I，计划(动态，16):0,000683s
• 并行 I，计划(静态，16):0,000647s
• 平行 J，无时间表:0,001978 秒(这是我执行的位置 执行速度较慢)

大小 = 256

• 序列号:0,005787s
• 并行 I，计划(动态，16):0,005125s
• 并行 I，计划(静态，16):0,004938s
• 并行 J，无调度:0,013916s

大小 = 1024

• 序列号:0,930250s
• 并行 I，计划(动态，16):0,865750s
• 并行 I，计划(静态，16):0,823750s
• 平行 J，无时间表:1,137000s

Answer 1

注意:此答案与如何从循环顺序中获得最佳性能或如何并行化它无关，因为出于多种原因我认为它不是最优的。我将尝试就如何改进顺序(并将其并行化)提供一些建议。

循环顺序

OpenMP 通常用于在多个 CPU 上分配工作。因此，您希望最大化每个线程的工作负载，同时最小化所需的数据和信息传输量。

1. 您想并行执行最外层循环而不是第二个循环。因此，您需要将 r_matrix 索引之一作为外循环索引，以避免在写入结果矩阵时出现竞争条件。

2. 接下来您要按内存存储顺序遍历矩阵(将变化较快的索引作为第二个而不是第一个下标索引)。

您可以通过以下循环/索引顺序实现两者:

for i = 0 to a_rows
  for k = 0 to a_cols
    for j = 0 to b_cols
      r[i][j] = a[i][k]*b[k][j]

在哪里

• j 比 i 或 k 变化快，k 比 i 变化快>.
• i是结果矩阵下标，i循环可以并行运行

以这种方式重新排列您的 multiply_matrices_KIJ 已经提供了相当多的性能提升。

我做了一些简短的测试，我用来比较时间的代码是:

template<class T>
void mm_kij(T const * const matrix_a, std::size_t const a_rows, 
  std::size_t const a_cols, T const * const matrix_b, std::size_t const b_rows, 
  std::size_t const b_cols, T * const matrix_r)
{
  for (std::size_t k = 0; k < a_cols; k++)
  {
    for (std::size_t i = 0; i < a_rows; i++)
    {
      for (std::size_t j = 0; j < b_cols; j++)
      {
        matrix_r[i*b_cols + j] += 
          matrix_a[i*a_cols + k] * matrix_b[k*b_cols + j];
      }
    }
  }
}

模仿您的multiply_matrices_KIJ() 函数对比

template<class T>
void mm_opt(T const * const a_matrix, std::size_t const a_rows, 
  std::size_t const a_cols, T const * const b_matrix, std::size_t const b_rows, 
  std::size_t const b_cols, T * const r_matrix)
{
  for (std::size_t i = 0; i < a_rows; ++i)
  { 
    T * const r_row_p = r_matrix + i*b_cols;
    for (std::size_t k = 0; k < a_cols; ++k)
    { 
      auto const a_val = a_matrix[i*a_cols + k];
      T const * const b_row_p = b_matrix + k * b_cols;
      for (std::size_t j = 0; j < b_cols; ++j)
      { 
        r_row_p[j] += a_val * b_row_p[j];
      }
    }
  }
}

执行上述命令。

Time consumption for multiplication of two 2048x2048 matrices on Intel i5-2500k

• mm_kij(): 6.16706s.

• mm_opt(): 2.6567s.

给定的顺序还允许外循环并行化，而不会在写入结果矩阵时引入任何竞争条件:

template<class T>
void mm_opt_par(T const * const a_matrix, std::size_t const a_rows, 
  std::size_t const a_cols, T const * const b_matrix, std::size_t const b_rows, 
  std::size_t const b_cols, T * const r_matrix)
{
#if defined(_OPENMP)
  #pragma omp parallel
  {
    auto ar = static_cast<std::ptrdiff_t>(a_rows);
    #pragma omp for schedule(static) nowait
    for (std::ptrdiff_t i = 0; i < ar; ++i)
#else
    for (std::size_t i = 0; i < a_rows; ++i)
#endif
    {
      T * const r_row_p = r_matrix + i*b_cols;
      for (std::size_t k = 0; k < b_rows; ++k)
      {
        auto const a_val = a_matrix[i*a_cols + k];
        T const * const b_row_p = b_matrix + k * b_cols;
        for (std::size_t j = 0; j < b_cols; ++j)
        {
          r_row_p[j] += a_val * b_row_p[j];
        }
      }
    }
#if defined(_OPENMP)
  }
#endif
}

每个线程写入一个单独的结果行

Time consumption for multiplication of two 2048x2048 matrices on Intel i5-2500k (4 OMP threads)

• mm_kij(): 6.16706s.

• mm_opt(): 2.6567s.

• mm_opt_par(): 0.968325s.

不是完美的缩放，但作为一个比串行代码更快的开始。