使用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]);
                }       
            }       
        }   
    }

您能告诉我我在这里缺少什么，或者为什么这不是一个劣等/优越的解决方案？

Answer 1

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

如果要编写可扩展的代码，则需要最小化同步和通信（包括原子操作）。结果，使用collapse(2)比collapse(3)好得多。但这不是唯一的问题是您的代码。实际上，要提高效率，您必须连续执行内存访问，并尽可能将数据保留在缓存中。

例如，由于更多连续的内存访问（在我的PC上大约是8倍，在大多数系统上，交换在i上迭代的循环和在k上迭代的循环（不适用于collapse(2)）要快几倍：] >

    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，OpenBLAS，ATLAS或Eigen之类的Intel MKL库，而不是编写自己的代码。实际上，此类库已经过非常优化，并且通常可以在许多内核上很好地扩展。如果您的目标是了解如何编写有效的矩阵乘法代码，则好的开始可能是阅读this tutorial。

Answer 2

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