用于矩阵乘法的 OpenMP答案

【问题标题】：OpenMP for matrix multiplication用于矩阵乘法的 OpenMP
【发布时间】：2014-08-15 11:26:25
【问题描述】：

我是 OpenMP 新手，正在拼命学习。我曾尝试在 Visual Studio 2012 中用 C++ 编写示例代码来实现矩阵乘法。我希望有 OpenMP 经验的人可以看看这段代码并帮助我获得最终的速度/并行化：

#include <iostream>
#include <stdlib.h>
#include <omp.h>
#include <random>
using namespace std;

#define NUM_THREADS 4

// Program Variables
double**        A;
double**        B;
double**        C;
double          t_Start;
double          t_Stop;
int             Am;
int             An;
int             Bm;
int             Bn;

// Program Functions
void            Get_Matrix();
void            Mat_Mult_Serial();
void            Mat_Mult_Parallel();
void            Delete_Matrix();


int main()
{
    printf("Matrix Multiplication Program\n\n");
    cout << "Enter Size of Matrix A: ";
    cin >> Am >> An;
    cout << "Enter Size of Matrix B: ";
    cin >> Bm >> Bn;

    Get_Matrix();
    Mat_Mult_Serial();
    Mat_Mult_Parallel();


    system("pause");
    return 0;

}


void Get_Matrix()
{
    A = new double*[Am];
    B = new double*[Bm];
    C = new double*[Am];
    for ( int i=0; i<Am; i++ ){A[i] = new double[An];}
    for ( int i=0; i<Bm; i++ ){B[i] = new double[Bn];}
    for ( int i=0; i<Am; i++ ){C[i] = new double[Bn]; }

    for ( int i=0; i<Am; i++ )
    {
         for ( int j=0; j<An; j++ )
         {
             A[i][j]= rand() % 10 + 1;
         }
    }

    for ( int i=0; i<Bm; i++ )
    {
        for ( int j=0; j<Bn; j++ )
        {
            B[i][j]= rand() % 10 + 1;
        }
    }
    printf("Matrix Create Complete.\n");
}


void Mat_Mult_Serial()
{
    t_Start = omp_get_wtime();
    for ( int i=0; i<Am; i++ )
    {
        for ( int j=0; j<Bn; j++ )
        {
            double temp = 0;
            for ( int k=0; k<An; k++ )
            {
                temp += A[i][k]*B[k][j];
            }
        }
    }
    t_Stop = omp_get_wtime() - t_Start;
    cout << "Serial Multiplication Time: " << t_Stop << " seconds" << endl;
    }


void Mat_Mult_Parallel()
{
    int i,j,k;
    t_Start = omp_get_wtime();

    omp_set_num_threads(NUM_THREADS);
    #pragma omp parallel for private(i,j,k) schedule(dynamic)
    for ( i=0; i<Am; i++ )
    {
        for ( j=0; j<Bn; j++ )
        {
            //double temp = 0;
            for ( k=0; k<An; k++ )
            {
                C[i][j] += A[i][k]*B[k][j];
            }
        }
    }

    t_Stop = omp_get_wtime() - t_Start;
    cout << "Parallel Multiplication Time: " << t_Stop << " seconds." << endl;
}


void Delete_Matrix()
{
    for ( int i=0; i<Am; i++ ){ delete [] A[i]; }
    for ( int i=0; i<Bm; i++ ){ delete [] B[i]; }
    for ( int i=0; i<Am; i++ ){ delete [] C[i]; }

    delete [] A;
    delete [] B;
    delete [] B;
}

【问题讨论】：

我有两个 cmets。首先是你可能不应该并行化k。由于您反复修改C[i][j]，我认为这些操作不能有效地并行化。（并行化i 和j 应该没问题）第二个是内存局部性和缓存未命中往往会在此类代码中产生最大的差异，因此您可能需要考虑存储B 的转置而不是@ 987654327@ 本身以获得最佳性能。（假设A 和B 很大）

标签： c++ matrix-multiplication openmp

【解决方案1】：

我对 OpenMP 很陌生，这段代码很有启发性。然而，我在串行版本中发现了一个错误，与并行版本相比，它具有不公平的速度优势。

您在串行版本中写了temp += A[i][k]*B[k][j];，而不是像在并行版本中那样写C[i][j] += A[i][k]*B[k][j];。这要快得多（但不能帮助您计算 C 矩阵）。所以你不是在比较苹果和苹果，这使得并行代码相比之下看起来更慢。当我修复这条线并在我的笔记本电脑上运行它（允许 2 个线程）时，并行版本的速度几乎快了两倍。还不错！

【讨论】：

【解决方案2】：

我的示例基于我为并行教学创建的矩阵类。如果您有兴趣，请随时与我联系。有几种方法可以加速矩阵乘法：

存储

按行主要顺序使用一维数组以更快地访问元素。
您可以使用 A[i * An + j] 访问 A(i,j)

使用循环不变优化

for (int i = 0; i < m; i ++)
    for (int j = 0; j < p; j ++)
    {
        Scalar sigma = C(i, j);
        for (int k = 0; k < n; k ++)
            sigma += (*this)(i, k) * B(k, j);
        C(i, j) = sigma;
    }

这可以防止在最内层循环中多次重新计算 C(i,j)。

改变循环顺序“for k for i”

for (int i = 0; i < m; i ++)
    for (int k = 0; k < n; k ++)
    {
        Aik = (*this)(i, k);
        for (int j = 0; j < p; j ++)
            C(i, j) += Aik * B(k, j);
    }

这允许使用空间data locality

使用循环阻塞/平铺

for(int ii = 0; ii < m; ii += block_size)
    for(int jj = 0; jj < p; jj += block_size)
        for(int kk = 0; kk < n; kk += block_size)
            #pragma omp parallel for // I think this is the best place for this case
            for(int i = ii; i < ii + block_size; i ++)
                for(int k = kk; k < kk + block_size; k ++)
                {
                    Scalar Aik = (*this)(i, k);
                    for(int j = jj; j < jj + block_size; j ++)
                        C(i, j) +=  Aik * B(k, j);
                }

这可以使用更好的时间数据局部性。最佳 block_size 取决于您的架构和矩阵大小。

然后并行化！

一般来说，#pragma omp parallel for 应该在最外层循环中完成。也许在两个第一个外部循环中使用两个并行循环可以提供更好的结果。这取决于您使用的架构、矩阵大小......您必须进行测试！由于矩阵乘法具有静态工作负载，我将使用静态调度。

摩尔优化！

你可以loop nest optimization。您可以矢量化您的代码。你可以看看BLAS是怎么做的。

【讨论】：