矩阵乘法、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 for 循环并检查性能提升。我的问题是，我期望（ed）至少快一点执行（比我现在得到的最多快 5-10%）即使它是并行的 I 循环（不能用K 循环，因为我会得到不正确的结果，因为它是 matrix_r[i][j])。

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

大小 = 128

• 系列版本：0,000608s
• 并行 I，调度（动态，16）：0,000683s
• 并行 I，调度（静态，16）：0,000647s
• 并行 J，无计划：0,001978s（这是我执行的地方 执行速度较慢）

大小 = 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];
      }
    }
  }
}

执行上述命令。

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

• 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
}

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

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

• mm_kij(): 6.16706s.

• mm_opt()：2.6567s。

• mm_opt_par(): 0.968325s.

不是完美的缩放，但比串行代码更快。

Answer 2

OpenMP 实现会创建一个线程池（尽管 OpenMP 标准并未强制要求线程池，但我见过的每个 OpenMP 实现都会这样做），这样就不必在每次进入并行区域时都创建和销毁线程.然而，每个并行区域之间存在障碍，因此所有线程都必须同步。在并行区域之间的分叉连接模型中可能存在一些额外的开销。因此，即使不必重新创建线程，它们仍然必须在并行区域之间进行初始化。更多详情请见here。

为了避免进入并行区域之间的开销，我建议在最外层循环上创建并行区域，但在i 上进行内部循环上的工作共享，如下所示：

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];
}

使用#pragma omp for 时存在隐式障碍。 nowait 子句消除了障碍。

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

Answer 3

请始终牢记，出于缓存目的，循环的最佳排序将是最慢的 -> 最快的。在您的情况下，这意味着 I,K,L 顺序。如果您的编译器没有从 KIJ->IKL 排序自动重新排序您的串行代码（假设您有“-O3”），我会感到非常惊讶。但是，编译器无法对您的并行循环执行此操作，因为这会破坏您在并行区域内声明的逻辑。

如果您真的无法重新排序循环，那么您最好的选择可能是重写并行区域以包含可能的最大循环。如果您有 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];
    }
}