矩阵乘法的平行优化
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
矩阵乘法的并行计算及其优化
郭善良∗ 上海师范大学数理学院 信息与计算科学专业
摘要
在本文中我们先将矩阵乘法进行并行化处理，然后再在此基础上对其进行了 几种优化，并对每种优化进行了算例的验证计算。
在数学研究和工程计算中，许多问题最后都可以归纳成为解一个常（偏）微分方程组，而 大部分此类方程都无法求得其解析解，所以通常是对其进行离散化处理后转化成一个线性系 统，通俗地讲就是解线性方程组或求线性迭代式的值。这其中便涉及到大量的矩阵运算，尤 其是矩阵的乘法运算。矩阵乘法看似简单但大矩阵的乘法却是非常消耗计算机的运算时间， 例如：一个3000阶方阵（这在工程上是非常常见的规模）做一次乘法一般要花时几十到上百秒 （看计算机的性能而定），如果矩阵的阶数再大上去其运算的增长量将是按矩阵阶的三次方关 系上升的，所以优化这个运算就非常有意义。 目前我们所使用的计算机不管是台式机还是手持机（平板和手机）都已经内置了好几个计 算核心，且通常我们的机器都带有近期的显卡，而这些显卡基本上都支持显卡计算。但我们平 时的编程还停留在单进（线）程的CPU计算阶段上，大量的现存资源被白白地浪费了。本文就 是讨论一下矩阵乘法的并行化处理及其优化问题。以使我们能充分利用手中的计算资源来提高 矩阵乘法的速度。 在本文中所有涉及到的计算的软件环境是：64位Windows7 SP1操作系统和Visual C++ 2013(Express)的编译器，编程语言为：C++、OpenMP和OpenCL；硬件环境是：Xeon X5672的 四核CPU（3.2G Hz，相当于第一代i7 CPU），18G内存，500G7200转的西部数据硬盘，Nvidia GTX580显卡（基本上是已经落后的硬件环境，这样配置的机器可以在二手市场上以非常便宜 的价格买到），编译的优化参数为-O2。
1 矩阵乘法的CPU并行化
假设我们有两个矩阵A = (aij )m×n 和B = (bij )n×k ，而C = (cij )m×k 为这两个矩阵的乘 积。则我们知道： cij =
∗ 学号:
n ∑ p=1
aip × bpj
xxxxxxxxxxxx
1
令数组a[m × n]，b[n × k ]，c[m × k ]为三个以行优先方式保存三个矩阵的数组。则上面的矩 阵乘法的伪代码为： const int nN=3000,nM=nN*nN; double *pA,*pB,*pC; //三个矩阵 ... for (int i = 0; i < nM; ++i) { int row = i / nN; int col = i%nN; double d = 0; for (int j = 0; j < nN; ++j) d += pA[row*nN + j] * pB[j*nN + col]; pC[i] = d; } 我们可以非常容易将上面的代码改写成OpenMP形式的并行伪代码为： #pragma omp parallel for firstprivate(pA,pB,pC) for (int i = 0; i < nM; ++i) { int row = i / nN; int col = i%nN; double d = 0; for (int j = 0; j < nN; ++j) d += pA[row*nN + j] * pB[j*nN + col]; pC[i] = d; } 经过算例验证，我们在矩阵为3000阶的时候运算时间从284.334202秒减少到65.712996秒。 这说明在多核计算机上即使通过非常简单的并行化处理就可使得矩阵乘法得到极大加速。
参考文献
[1] Peter S. Pacheco，并行程序设计导论，机械工业出版社，2012，8 [2] 周伟民，OpenMP简易教程，百度文库 [3] Gene H. Golub，矩阵计算（第三版），人民邮电出版社，2012，6 [4] Nvidia，CUBLAS Library，Nvidia，2014,8 [5] Nvidia，CUDA C Programming Guide，Nvidia，2014,8
3 结语
由于矩阵的乘法的应用的广泛性，所有只要运算效率有所提高都是有意义的，而本文中的 优化使得在低端的双4核（8核）计算机上且与免费的计算软件配合下，其运算效率提高了26倍 以上，而应用GTX580显卡来加速更可以将效率提高到360倍以上，如果计算机中CPU的计算核 心更多或配置一块更新更好的显卡或专业计算卡，其效率将还要高，所以我们对此是感到非常 地满意，而且特别是在CPU的并行优化中对第二个矩阵进行转置的优化是本文第一次提出的， 有一定的创新性，应用这种优化哪怕不进行并行处理都可以使计算效率大大地提高。
1.1 算法的优化
我们知道，现在的CPU都在内部集成了高速缓存，所以访问内存中连续的数据时，速度要 远快于不连续的数据。而上面代码中对矩阵A的数据访问是连续的，但对矩阵B 的数据访问却是 不连续，所以我们在做乘法前先对矩阵B 作转置处理，这样对B 的访问也是连续的。具体修改代 码如下： double *pD;//用作临时矩阵 ... #pragma omp parallel firstprivate(pA,pB,pC)
} double *pS=(double*)&psseS; dT=pS[0]+pS[1]; if((n & 1)!=0)//计算后面遗留的一个 { p1=(double *)psse_a; p2=(double *)psse_d; dT+=(*p1)*(*p2); } pC[k]=dT; } } 我们再通过具体算例验算得，运算的时间为：4.481360秒，速度有了进一步的提高，但编 程的难度也有提高，当然我们还有一条优化路径未试，即通过汇编指令来使用多媒体指令集， 但这个编程难度更是高许多，且提升的速度估计不会比用C++方式使用SSE2指令集快许多。 这样我们把需花时284.334202秒的3000阶矩阵乘法运算优化到4.481360秒，提高了计算 效率63倍以上，接着我们试算了一下原来似乎等不到尽头的5000阶矩阵的乘法，现在大概 在38.607731秒左右就可以得到结果，这使得我们更有信心来计算大数据的矩阵乘法了。另 外，如果资金充足的话还可以向Intel购买Intel C++以及Intel MKL（数学并行运算库，但这个 软件价格不便宜啊），应用其中的cblas库函数来做矩阵的乘法可以使得运算的时间更为缩短 （5000阶矩阵的乘法大概在5.286354秒左右完成）。
5
cublasGetMatrix(nN,nN, sizeof(double), d_pfC, nN, pfC, nN); 我们通过具体算例验算得，运算的时间为：0.34264秒，速度有10多倍以上的提高，从此处 我们看出显卡的运算能力要比CPU的运算强大的多。最后我们再次试算了一下5000阶矩阵的乘 法，现在大概在1.48494秒左右就可以得到结果。但是由于我们机器上的显卡是为游戏而设计， 不是专业的计算卡，所以它同时承担屏幕的刷新工作。如果你的显卡计算程序运行时间过长会 造成出错，所以如果我们的程序遇到这种情况，我们可以把大矩阵分块多次递交运算任务以便 让显卡有机会进行屏幕刷新工作。
{ d += (*p1) * (*p2); ++p1; ++p2; } pC[i] = d; } } 我们再通过具体算例验算得，运算的时间为：5.646863秒，优化程度不大，所以这样的优 化可以忽略。 现在的CPU都带有SSE多媒体快速运算指令集，所以我们可以用CPU的单指令多操作数 的SSE2多媒体指令集继续对矩阵乘法进行优化，这也是一种并行处理，即并行上加并行。具体 代码如下： #pragma omp parallel shared(pA,pB,pC,pD) { #pragma omp for for (int i=0;i<m*n;++i)//矩阵转置 { int r=i/n; int c=i%n; pD[c*n+r]=pB[i]; } #pragma omp for for (int i=0;i<m*k;++i) { int r=i/k; int c=i%k; double dT=0.0; double *p1=pA+r*n,*p2=pD+c*n; __m128d *psse_a=(__m128d*)p1; __m128d *psse_d=(__m128d*)p2; __m128d psseS=_mm_setzero_pd(); //计算前面成对的乘法，一次算二对双精度浮点数的乘法 for(int l=0;l<n/2;++l) { __m128d m=_mm_mul_pd(*psse_a,*psse_d); psseS=_mm_add_pd(psseS,m); ++psse_a; ++psse_d; 4
2
{ #pragma omp for for (int i = 0; i < nM; ++i)//B矩阵转置 { int row = i / nN; int col = i%nN; pD[col*nN + row] = pB[i]; } #pragma omp for for (int i = 0; i < nM; ++i) { int row = i / nN; int col = i%nN; double d = 0; for (int j = 0; j < nN; ++j) d += pA[row*nN + j] * pD[col*nN + j]; pC[i] = d; } } 通过具体算例验算得，此时加上矩阵转置一起3000阶矩阵的乘法所花时间比前面最快的代 码运行时间还要快，花时：5.765069秒。 在上面的代码中对矩阵的元素的访问都是通过pA[r*n+j]和pB[c*k+j]来实现的，在这里我 们反复进行了乘法的计算，所以我们还可以优化一下，优化的代码如下： #pragma omp parallel shared(m,n,k,pA,pB,pC,pD) { #pragma omp for for (int i = 0; i < nM; ++i)//B矩阵转置 { int row = i / nN; int col = i%nN; pD[col*nN + row] = pB[i]; } #pragma omp for for (int i = 0; i < nM; ++i) { int row = i / nN; int col = i%nN; double d = 0，*p1=pA+row*nN,*p2=pD+col*nN; for (int j = 0; j < nN; ++j) 3
2 矩阵乘法的GPU并行化
目前最新的高端显卡都具有上千个运算单元，所以单张显卡的浮点计算能力要远远大 于CPU的浮点运算能力，所以对矩阵乘法的并行优化除了可以进行CPU的并行优化外还可 以应用计算机的显卡上那数量恐怖的浮点运算单元来优化并行计算。应用Nvidia公司所提供 的CUDA库的运算代码如下： //在显卡上分配空间 cudaMalloc((void**)&d_pfA, nN*nN*sizeof(double)); cudaMalloc((void**源自文库&d_pfB, nN*nN*sizeof(double)); cudaMalloc((void**)&d_pfC, nN*nN*sizeof(double)); //将数据传送到显卡上 cublasSetMatrix(nN, nN, sizeof(double), pfA, nN, d_pfA, nN); cublasSetMatrix(nN, nN, sizeof(double), pfB, nN, d_pfB, nN); //计算 cublasDgemm(cubHandle, cubTrans,cubTrans,nN, nN,nN, &fAlpha, d_pfA, nN, d_pfB, nN, &fBeta, d_pfC, nN); //将数据传送回主机内存
郭善良∗ 上海师范大学数理学院 信息与计算科学专业
摘要
在本文中我们先将矩阵乘法进行并行化处理，然后再在此基础上对其进行了 几种优化，并对每种优化进行了算例的验证计算。
在数学研究和工程计算中，许多问题最后都可以归纳成为解一个常（偏）微分方程组，而 大部分此类方程都无法求得其解析解，所以通常是对其进行离散化处理后转化成一个线性系 统，通俗地讲就是解线性方程组或求线性迭代式的值。这其中便涉及到大量的矩阵运算，尤 其是矩阵的乘法运算。矩阵乘法看似简单但大矩阵的乘法却是非常消耗计算机的运算时间， 例如：一个3000阶方阵（这在工程上是非常常见的规模）做一次乘法一般要花时几十到上百秒 （看计算机的性能而定），如果矩阵的阶数再大上去其运算的增长量将是按矩阵阶的三次方关 系上升的，所以优化这个运算就非常有意义。 目前我们所使用的计算机不管是台式机还是手持机（平板和手机）都已经内置了好几个计 算核心，且通常我们的机器都带有近期的显卡，而这些显卡基本上都支持显卡计算。但我们平 时的编程还停留在单进（线）程的CPU计算阶段上，大量的现存资源被白白地浪费了。本文就 是讨论一下矩阵乘法的并行化处理及其优化问题。以使我们能充分利用手中的计算资源来提高 矩阵乘法的速度。 在本文中所有涉及到的计算的软件环境是：64位Windows7 SP1操作系统和Visual C++ 2013(Express)的编译器，编程语言为：C++、OpenMP和OpenCL；硬件环境是：Xeon X5672的 四核CPU（3.2G Hz，相当于第一代i7 CPU），18G内存，500G7200转的西部数据硬盘，Nvidia GTX580显卡（基本上是已经落后的硬件环境，这样配置的机器可以在二手市场上以非常便宜 的价格买到），编译的优化参数为-O2。
1 矩阵乘法的CPU并行化
假设我们有两个矩阵A = (aij )m×n 和B = (bij )n×k ，而C = (cij )m×k 为这两个矩阵的乘 积。则我们知道： cij =
∗ 学号:
n ∑ p=1
aip × bpj
xxxxxxxxxxxx
1
令数组a[m × n]，b[n × k ]，c[m × k ]为三个以行优先方式保存三个矩阵的数组。则上面的矩 阵乘法的伪代码为： const int nN=3000,nM=nN*nN; double *pA,*pB,*pC; //三个矩阵 ... for (int i = 0; i < nM; ++i) { int row = i / nN; int col = i%nN; double d = 0; for (int j = 0; j < nN; ++j) d += pA[row*nN + j] * pB[j*nN + col]; pC[i] = d; } 我们可以非常容易将上面的代码改写成OpenMP形式的并行伪代码为： #pragma omp parallel for firstprivate(pA,pB,pC) for (int i = 0; i < nM; ++i) { int row = i / nN; int col = i%nN; double d = 0; for (int j = 0; j < nN; ++j) d += pA[row*nN + j] * pB[j*nN + col]; pC[i] = d; } 经过算例验证，我们在矩阵为3000阶的时候运算时间从284.334202秒减少到65.712996秒。 这说明在多核计算机上即使通过非常简单的并行化处理就可使得矩阵乘法得到极大加速。
参考文献
[1] Peter S. Pacheco，并行程序设计导论，机械工业出版社，2012，8 [2] 周伟民，OpenMP简易教程，百度文库 [3] Gene H. Golub，矩阵计算（第三版），人民邮电出版社，2012，6 [4] Nvidia，CUBLAS Library，Nvidia，2014,8 [5] Nvidia，CUDA C Programming Guide，Nvidia，2014,8
3 结语
由于矩阵的乘法的应用的广泛性，所有只要运算效率有所提高都是有意义的，而本文中的 优化使得在低端的双4核（8核）计算机上且与免费的计算软件配合下，其运算效率提高了26倍 以上，而应用GTX580显卡来加速更可以将效率提高到360倍以上，如果计算机中CPU的计算核 心更多或配置一块更新更好的显卡或专业计算卡，其效率将还要高，所以我们对此是感到非常 地满意，而且特别是在CPU的并行优化中对第二个矩阵进行转置的优化是本文第一次提出的， 有一定的创新性，应用这种优化哪怕不进行并行处理都可以使计算效率大大地提高。
1.1 算法的优化
我们知道，现在的CPU都在内部集成了高速缓存，所以访问内存中连续的数据时，速度要 远快于不连续的数据。而上面代码中对矩阵A的数据访问是连续的，但对矩阵B 的数据访问却是 不连续，所以我们在做乘法前先对矩阵B 作转置处理，这样对B 的访问也是连续的。具体修改代 码如下： double *pD;//用作临时矩阵 ... #pragma omp parallel firstprivate(pA,pB,pC)
} double *pS=(double*)&psseS; dT=pS[0]+pS[1]; if((n & 1)!=0)//计算后面遗留的一个 { p1=(double *)psse_a; p2=(double *)psse_d; dT+=(*p1)*(*p2); } pC[k]=dT; } } 我们再通过具体算例验算得，运算的时间为：4.481360秒，速度有了进一步的提高，但编 程的难度也有提高，当然我们还有一条优化路径未试，即通过汇编指令来使用多媒体指令集， 但这个编程难度更是高许多，且提升的速度估计不会比用C++方式使用SSE2指令集快许多。 这样我们把需花时284.334202秒的3000阶矩阵乘法运算优化到4.481360秒，提高了计算 效率63倍以上，接着我们试算了一下原来似乎等不到尽头的5000阶矩阵的乘法，现在大概 在38.607731秒左右就可以得到结果，这使得我们更有信心来计算大数据的矩阵乘法了。另 外，如果资金充足的话还可以向Intel购买Intel C++以及Intel MKL（数学并行运算库，但这个 软件价格不便宜啊），应用其中的cblas库函数来做矩阵的乘法可以使得运算的时间更为缩短 （5000阶矩阵的乘法大概在5.286354秒左右完成）。
5
cublasGetMatrix(nN,nN, sizeof(double), d_pfC, nN, pfC, nN); 我们通过具体算例验算得，运算的时间为：0.34264秒，速度有10多倍以上的提高，从此处 我们看出显卡的运算能力要比CPU的运算强大的多。最后我们再次试算了一下5000阶矩阵的乘 法，现在大概在1.48494秒左右就可以得到结果。但是由于我们机器上的显卡是为游戏而设计， 不是专业的计算卡，所以它同时承担屏幕的刷新工作。如果你的显卡计算程序运行时间过长会 造成出错，所以如果我们的程序遇到这种情况，我们可以把大矩阵分块多次递交运算任务以便 让显卡有机会进行屏幕刷新工作。
{ d += (*p1) * (*p2); ++p1; ++p2; } pC[i] = d; } } 我们再通过具体算例验算得，运算的时间为：5.646863秒，优化程度不大，所以这样的优 化可以忽略。 现在的CPU都带有SSE多媒体快速运算指令集，所以我们可以用CPU的单指令多操作数 的SSE2多媒体指令集继续对矩阵乘法进行优化，这也是一种并行处理，即并行上加并行。具体 代码如下： #pragma omp parallel shared(pA,pB,pC,pD) { #pragma omp for for (int i=0;i<m*n;++i)//矩阵转置 { int r=i/n; int c=i%n; pD[c*n+r]=pB[i]; } #pragma omp for for (int i=0;i<m*k;++i) { int r=i/k; int c=i%k; double dT=0.0; double *p1=pA+r*n,*p2=pD+c*n; __m128d *psse_a=(__m128d*)p1; __m128d *psse_d=(__m128d*)p2; __m128d psseS=_mm_setzero_pd(); //计算前面成对的乘法，一次算二对双精度浮点数的乘法 for(int l=0;l<n/2;++l) { __m128d m=_mm_mul_pd(*psse_a,*psse_d); psseS=_mm_add_pd(psseS,m); ++psse_a; ++psse_d; 4
2
{ #pragma omp for for (int i = 0; i < nM; ++i)//B矩阵转置 { int row = i / nN; int col = i%nN; pD[col*nN + row] = pB[i]; } #pragma omp for for (int i = 0; i < nM; ++i) { int row = i / nN; int col = i%nN; double d = 0; for (int j = 0; j < nN; ++j) d += pA[row*nN + j] * pD[col*nN + j]; pC[i] = d; } } 通过具体算例验算得，此时加上矩阵转置一起3000阶矩阵的乘法所花时间比前面最快的代 码运行时间还要快，花时：5.765069秒。 在上面的代码中对矩阵的元素的访问都是通过pA[r*n+j]和pB[c*k+j]来实现的，在这里我 们反复进行了乘法的计算，所以我们还可以优化一下，优化的代码如下： #pragma omp parallel shared(m,n,k,pA,pB,pC,pD) { #pragma omp for for (int i = 0; i < nM; ++i)//B矩阵转置 { int row = i / nN; int col = i%nN; pD[col*nN + row] = pB[i]; } #pragma omp for for (int i = 0; i < nM; ++i) { int row = i / nN; int col = i%nN; double d = 0，*p1=pA+row*nN,*p2=pD+col*nN; for (int j = 0; j < nN; ++j) 3
2 矩阵乘法的GPU并行化
目前最新的高端显卡都具有上千个运算单元，所以单张显卡的浮点计算能力要远远大 于CPU的浮点运算能力，所以对矩阵乘法的并行优化除了可以进行CPU的并行优化外还可 以应用计算机的显卡上那数量恐怖的浮点运算单元来优化并行计算。应用Nvidia公司所提供 的CUDA库的运算代码如下： //在显卡上分配空间 cudaMalloc((void**)&d_pfA, nN*nN*sizeof(double)); cudaMalloc((void**源自文库&d_pfB, nN*nN*sizeof(double)); cudaMalloc((void**)&d_pfC, nN*nN*sizeof(double)); //将数据传送到显卡上 cublasSetMatrix(nN, nN, sizeof(double), pfA, nN, d_pfA, nN); cublasSetMatrix(nN, nN, sizeof(double), pfB, nN, d_pfB, nN); //计算 cublasDgemm(cubHandle, cubTrans,cubTrans,nN, nN,nN, &fAlpha, d_pfA, nN, d_pfB, nN, &fBeta, d_pfC, nN); //将数据传送回主机内存