需要帮助了解OpenMP矩阵乘法C＆＃x2B;＆＃x2B;代码

Question

这是我编写的矩阵乘法C ++ OpenMP代码。我正在尝试使用OpenMP来优化程序。顺序代码速度为7秒，但是当我添加OpenMP语句时，但是它的速度仅为3秒。我认为它会变得更快，不明白我是否做对了。

OpenMP语句位于fill_random函数中，并且在MAIN中的循环部分的矩阵乘法三重。

我将感谢您可以提供的任何帮助或建议，以理解这一点！

#include <iostream>
#include <cassert>
#include <omp.h>
#include <chrono>

using namespace std::chrono;


double** fill_random(int rows, int cols )
{
    
    double** mat = new double* [rows]; //Allocate rows.
    #pragma omp  parallell collapse(2) 
    for (int i = 0; i < rows; ++i)
    {
        mat[i] = new double[cols];           // added
        for( int j = 0;  j < cols; ++j)
        {
            mat[i][j] = rand() % 10;
        }
       
    }
     return mat;
}


double** create_matrix(int rows, int cols)
{
    double** mat = new double* [rows]; //Allocate rows.
    for (int i = 0; i < rows; ++i)
    {
        mat[i] = new double[cols](); //Allocate each row and zero initialize..
    }
    return mat;
}

void destroy_matrix(double** &mat, int rows)
{
    if (mat)
    {
        for (int i = 0; i < rows; ++i)
        {
            delete[] mat[i]; //delete each row..
        }

        delete[] mat;  //delete the rows..
        mat = nullptr;
    }
}

int main()
{
    int rowsA = 1000; // number of rows
    int colsA= 1000; // number of columns
    double** matA = fill_random(rowsA, colsA);


    int rowsB = 1000; // number of rows
    int colsB = 1000; // number of columns
    double** matB = fill_random(rowsB, colsB);


//Checking matrix multiplication qualification
    assert(colsA == rowsB);


    double** matC = create_matrix(rowsA, colsB);

    //measure the multiply only
    const auto start = high_resolution_clock::now();

    //Multiplication
    #pragma omp parallel for 
    
    for(int i = 0; i < rowsA; ++i)
    {
        for(int j = 0; j < colsB; ++j)
        {
            for(int k = 0; k < colsA; ++k) //ColsA..
            {
                matC[i][j] += matA[i][k] * matB[k][j];
            }
        }
        
    }

    const auto stop = high_resolution_clock::now();
    const auto duration = duration_cast<seconds>(stop - start);

    std::cout << "Time taken by function: " << duration.count() << " seconds" << std::endl;



    //Clean up..
    destroy_matrix(matA, rowsA);
    destroy_matrix(matB, rowsB);
    destroy_matrix(matC, rowsA);

    return 0;
}

Answer 1

1. 您的问题很小。
2. 矩阵创建中的崩溃，因为循环不是完美嵌套的，因此无能为力。另一方面，在乘法例程中，您应添加collapse（2）指令。
3. 用一系列指针创建一个矩阵意味着表达式matb [k] [j]在内存中跳舞。将矩阵分配为一个数组，然后使用i*n+j作为索引表达式。 （当然，我会把它放在宏观上。）

Answer 2

具有双（64位）元素类型的1000x1000的矩阵大小需要8MB数据。当您乘以两个矩阵时，您会读取16MB数据。当您写入第三个矩阵时，您还可以访问24MB数据。

如果L3缓存小于24MB，则RAM是瓶颈。也许单线线程没有完全使用其带宽，但是当使用OpenMP时，RAM带宽将完全使用。就您而言，带宽只有50％的余量。

幼稚版本的使用良好。您需要交换两个循环的顺序以获得更多的缓存：

loop
   loop k
     loop
       C[..] += B[..] * A[..]

尽管增加C不使用此优化版本中的寄存器，但在这种情况下，它重新使用缓存，这在这种情况下更为重要。如果这样做，即使在单线程中，它也应该获得约100-200毫秒的计算时间。

另外，如果您需要性能，请不要这样做：

//Allocate each row and zero initialize..

一次分配整个矩阵，以使您的矩阵不会散布在内存中。

为了有效地添加更多线程，您可以执行子矩阵乘法以计算完整的矩阵乘法。扫描线乘法不利于线程之间的负载平衡。当子矩阵乘以时，由于缓存和从内存中获取的每个元素的浮点操作更高的浮点操作，它们会提供更好的负载分布。

编辑：

循环的交换顺序还使编译器能够矢量化最终的循环，因为其中一个输入矩阵在最内向的循环中变为常数。