矩阵乘法，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社区版，在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，并检查性能的提高。我的问题是我期望至少有一点点更快的执行速度(比我现在得到的最多快5-10%)，即使它是并行的I循环(不能用K循环，因为它是matrix_r[i][j])。

这些是我在使用上面显示的代码时得到的结果（我进行了数百次计算，得到了平均时间）：

• 序列版：0,000608s
• 并行I，时间表（动态，16）：0,000683s
• 并行I，时间表（静态，16）：0,000647s
• 并行J，没有计划：0,001978s（这是我遇到执行速度较慢的地方）

• 串行版本： 0，005787s
• 并行 I，计划（动态，16）：0，005125s
• 并行 I，时间表（静态，16）：0，004938s
• 并行J，无时间表：0，013916s

• 串行版本： 0，930250s
• 并行 I，时间表（动态，16）：0，865750s
• 并行 I，时间表（静态，16）：0，823750s
• 并行J，无时间表：1，137000s

Answer 1

永远记住，出于缓存目的，循环的最佳排序将是最慢的 -

如果你真的不能重新排序你的循环，那么你最好的选择可能是重写并行区域，以包含尽可能大的循环。如果你有OpenMP 4.0，你也可以考虑在你最快的维度上使用SIMD矢量化。然而，我仍然怀疑你是否能够击败你的串行代码，因为你的KIJ排序中固有的上述缓存问题…

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

Answer 2

OpenMP实现创建了一个线程池(尽管OpenMP标准没有强制要求线程池，但我见过的每个OpenMP实现都是这样做的),这样就不必在每次进入并行区域时都创建和销毁线程。然而，每个并行区域之间有一个屏障，因此所有线程必须同步。在并行区域之间的fork join模型中可能会有一些额外的开销。因此，即使不需要重新创建线程，它们仍然需要在并行区域之间进行初始化。更多细节可以在这里找到。

为了避免进入并行区域之间的开销，我建议在最外层的循环上创建并行区域，但在内部循环上进行工作共享，我喜欢这样：

void multiply_matrices_KIJ() {
    #pragma omp parallel
    for (int k = 0; k < SIZE; k++)
        #pragma omp for schedule(static) 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];
}

使用< code>#pragma omp for时存在一个隐含的障碍。< code>nowait子句移除该障碍。

还要确保您使用优化进行编译。在未启用优化的情况下比较性能几乎没有意义。我会使用-O3。

Answer 3

注意：这个答案不是关于如何从循环顺序中获得最佳性能或如何并行化它，因为我认为由于多种原因，它是次优的。我将尝试就如何改进顺序（并并行化它）提供一些建议。

循环顺序

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

>

您希望并行执行最外层的循环，而不是第二个循环。因此，您需要将其中一个r_matrix索引作为外部循环索引，以避免在写入结果矩阵时出现争用条件。

接下来的事情是，您希望按照内存存储顺序遍历矩阵（将变化更快的索引作为第二个而不是第一个下标索引）。

您可以使用以下循环/索引顺序实现这两个目标：

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]

在哪里

< li> j比< code>i或< code>k变化快，< code>k比< code>i变化快。 < li> i是结果矩阵下标，并且< code>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];
      }
    }
  }
}

执行上述命令。

Intel i5-2500k上两个2048x2048矩阵相乘的时间消耗

>

mm_kij（）： 6.16706s。

mm_opt（）：2.6567秒。

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

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
}

每个线程写入单个结果行的位置

在英特尔 i5-2500k 上乘法两个 2048x2048 矩阵（4 个 OMP 线程）所消耗的时间

>

mm_kij（）： 6.16706s。

mm_opt（）：2.6567秒。

< code>mm_opt_par()，0.968325s

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