使用 OpenMP 的矩阵乘法 (C) - 折叠所有循环

Question

所以我学习了 C 中的 OpenMP 基础知识和工作共享构造，尤其是循环。所有教程中使用的最著名示例之一是矩阵乘法，但它们都只是并行化了一个外循环或两个外循环。我想知道为什么我们不像我在这里所做的那样并行化和折叠所有 3 个循环（使用原子）：

    for(int i=0;i<100;i++){
        //Initialize the arrays
        for(int j=0;j<100;j++){
            A[i][j] = i;
            B[i][j] = j;
            C[i][j] = 0;       
        }       
    }

    //Starting the matrix multiplication
    #pragma omp parallel num_threads(4)
    {
        #pragma omp for collapse(3)
        for(int i=0;i<100;i++){
            for(int j=0;j<100;j++){
                for(int k=0;k<100;k++){
                        #pragma omp atomic
                        C[i][j] = C[i][j]+ (A[i][k]*B[k][j]);
                }       
            }       
        }   
    }

你能告诉我我在这里遗漏了什么或者为什么这不是 inferior/superior 解决方案吗？

Answer 1

与非原子操作相比，原子操作在大多数体系结构上的成本非常高（有关更详细的分析，请参阅 here to understand why or here）。当许多线程对同一共享内存区域进行并发访问时尤其如此。简而言之，一个原因是由于来自缓存一致性协议的隐式同步和通信，执行原子操作的线程不能在大多数硬件上不等待其他线程的情况下完全运行并行。减速的另一个来源是原子操作的高延迟（同样是由于缓存层次结构）。

如果您想编写可伸缩的代码，您需要尽量减少同步和通信（包括原子操作）。因此，使用 collapse(2) 比 collapse(3) 好得多。但这不是您的代码的唯一问题。事实上，为了提高效率，您必须连续执行内存访问并将数据尽可能多地保存在缓存中。

例如，交换遍历 i 的循环和遍历 k 的循环（不适用于 collapse(2)）在大多数系统上由于更连续的内存访问（大约 8 倍）在大多数系统上快几倍我的电脑）：

    for(int i=0;i<100;i++){
        //Initialize the arrays
        for(int j=0;j<100;j++){
            A[i][j] = i;
            B[i][j] = j;
            C[i][j] = 0;       
        }       
    }

    //Starting the matrix multiplication
    #pragma omp parallel num_threads(4)
    {
        #pragma omp for
        for(int i=0;i<100;i++){
            for(int k=0;k<100;k++){
                for(int j=0;j<100;j++){
                    C[i][j] = C[i][j] + (A[i][k]*B[k][j]);
                }
            }
        }
    }

编写快速矩阵乘法代码并不容易。如果您的目标是在生产代码中使用它，请考虑使用 BLAS libraries such as OpenBLAS, ATLAS, Eigen or Intel MKL 而不是编写您自己的代码。事实上，这样的库经过了非常优化，并且通常可以在许多内核上很好地扩展。如果您的目标是了解如何编写高效的矩阵乘法代码，那么阅读 this tutorial.

可能是一个很好的起点

Answer 2

折叠循环要求您知道自己在做什么，因为它可能会导致对缓存非常不友好的迭代拆分 space 或引入数据依赖性，具体取决于循环计数的乘积与数字的关系线程数。

想象一下下面构造的示例，这实际上并不少见（循环次数很少只是为了说明这一点）：

for (int i = 0; i < 7; i++)
  for (int j = 0; j < 3; j++)
     a[i] += b[i][j];

如果将外循环并行化，三个线程会进行两次迭代，一个线程只会进行一次迭代，但是所有线程都会执行内循环的所有迭代：

---0-- ---1-- ---2-- -3- (thread number)
000111 222333 444555 666 (values of i)
012012 012012 012012 012 (values of j)

每个 a[i] 仅由一个线程处理。聪明的编译器可能会使用寄存器优化来实现内部循环，首先将值累加到寄存器中，最后才分配给a[i]，它会非常快运行。

如果你打破这两个循环，你最终会遇到一个非常不同的情况。由于现在总共有 7x3 = 21 次迭代，默认拆分将是（取决于编译器和 OpenMP 运行时间，但大多数都是这样做的）每个线程五次迭代，一个线程进行六次迭代：

--0-- --1-- --2-- ---3-- (thread number)
00011 12223 33444 555666 (values of i)
01201 20120 12012 012012 (values of j)

如您所见，现在 a[1] 由线程 0 和线程 1 处理。类似地，a[3] 由线程 1 和线程 2 处理。这就是它 - 你刚刚引入了前一个案例中不存在的数据依赖性，因此现在您必须使用 atomic 来防止数据竞争。您为同步付出的代价比或多或少地进行一次迭代要高得多！在您的情况下，如果您只折叠两个外部循环，则根本不需要使用 atomic （尽管在您的特定情况下，4 除以 100，即使将所有循环折叠在一起，您也不需要需要 atomic 构造，但在一般情况下你需要它）。

另一个问题是，在折叠循环后，只有一个循环索引，i 和 j 索引都必须使用除法和模运算从这个新索引中重建。对于像您这样的简单循环体，重建索引的开销可能太高了。

Answer 3

没有什么好的理由不使用库进行矩阵-矩阵乘法，所以正如已经建议的那样，请调用 BLAS 而不是自己编写。尽管如此，您提出的问题并不特定于矩阵-矩阵乘法，因此无论如何都值得回答。

这里有几点可以改进：

使用连续内存。
如果 K 是最内层的循环，那么您正在做点积，这更难向量化。例如，循环顺序 IKJ 将更好地向量化。
如果您想使用 OpenMP 并行化点积，请使用归约而不是许多原子。

我在下面分别说明了这些技术。

连续内存

int n = 100;
double * C = malloc(n*n*sizeof(double));
for(int i=0;i<n;i++){
  for(int j=0;j<n;j++){
    C[i*n+j] = 0.0;
  }       
}

IKJ 循环顺序

for(int i=0;i<100;i++){
  for(int k=0;k<100;k++){
    for(int j=0;j<100;j++){
      C[i][j] = C[i][j]+ (A[i][k]*B[k][j]);
    }
  }
}

平行点积

double x = 0;
#pragma omp parallel for reduction(+:x)
for(int k=0;k<100;k++){
  x += (A[i][k]*B[k][j]);
}
C[i][j] += x;

外部资源

如何快速编写数字代码：简短介绍更详细地介绍了这些主题。

BLISlab 是专门针对矩阵-矩阵乘法的优秀教程，它将教您专家如何编写 BLAS 库调用。

使用 OpenMP 的矩阵乘法 (C) - 折叠所有循环

Matrix Multiplication using OpenMP (C) - Collapsing all the loops

c

parallel-processing

performance

openmp

matrix-multiplication

连续内存

IKJ 循环顺序

平行点积

外部资源