矩阵乘法的并行计算

转载：http://blog.csdn.net/realxie/article/details/7260072

设两个矩阵A和B，大小分别为M * N 和 N * P， 如果C = A * B， 则C的大小为M * P。

矩阵算法的算法表示，伪代码如下：

1. for (i = 0; i < M; ++i){  
2.     for (j = 0; j < P; ++j){  
3.         C[i][j] = 0;  
4.         for (k = 0; k < N; ++k){  
5.             C[i][j] += A[i][k] * B[k][j];  
6.         }  
7.     }  
8. }  

从上面的算法中可以看出该算法的时间复杂度为O(M*N*P)，当M,N,P都非常大时该计算将非常耗时。那么如何将上面的串行算法转换成并行算法呢？

从上面的三层循环中可以看出最外层的循环是独立的，即对C[i][*]的计算不依赖于任何C[ii][*]的计算，因此我们可以非常容易将最外层的循环转换成并行。

1. #prama omp parallel for num_threads(CORE_NUM)  
2. for (i = 0; i < M; ++i){  
3.     for (j = 0; j < P; ++j){  
4.         C[i][j] = 0;  
5.         for (k = 0; k < N; ++k){  
6.             C[i][j] += A[i][k] * B[k][j];  
7.         }  
8.     }  
9. }  

但是这里有一个局限，如果假设cpu的核数CORE_NUM > M，同样无法充分利用所有的计算资源。

进一步分析， 由于C矩阵的大小为M * P，那么我们能不能将C的计算下平均分配到CORE_NUM个核心上呢，即每个核分配ceil(M*P/CORE_NUM)个计算任何，即将上面的第一和第二层并行化。

首先将C转换成一维的数组T[M*P] , 则C[i][j] = T[i * M + j], 反过来T[z] = C[z/M] [ z %P]。

故进一步的并行算法为：

1. #prama omp parallel for num_threads(NUM)  
2. for (z = 0; z < M * P; ++z){  
3.         i = z / P;  
4.         j = z % P;  
5.         C[i][j] = 0;  
6.         for (k = 0; k < N; ++k){  
7.             C[i][j] += A[i][k] * B[k][j];  
8.     }  
9. }  

性能优化。

看最里面一层的计算

1. for (k = 0; k < N; ++k){  
2.     C[i][j] += A[i][k] * B[k][j];  

由于内存中二维数组是以行优先进行存储的，因此B[k][j]存在严重的cache命中率问题，解决这个问题的方法是也将B进行一次沿对角线进行翻转，使得最里面的计算变成

1. for (k = 0; k < N; ++k){  
2.     C[i][j] += A[i][k] * B[j][k];  

另外一点需要注意的就是C[i][j] += A[i][k] * B[j][k];计算时的伪共享问题。