设两个矩阵A和B,大小分别为M * N 和 N * P, 如果C = A * B, 则C的大小为M * P。
矩阵算法的算法表示,伪代码如下:
for (i = 0; i < M; ++i){
for (j = 0; j < P; ++j){
C[i][j] = 0;
for (k = 0; k < N; ++k){
C[i][j] += A[i][k] * B[k][j];
}
}
}
从上面的算法中可以看出该算法的时间复杂度为O(M*N*P),当M,N,P都非常大时该计算将非常耗时。那么如何将上面的串行算法转换成并行算法呢?
从上面的三层循环中可以看出最外层的循环是独立的,即对C[i][*]的计算不依赖于任何C[ii][*]的计算,因此我们可以非常容易将最外层的循环转换成并行。
#prama omp parallel for num_threads(CORE_NUM)
for (i = 0; i < M; ++i){
for (j = 0; j < P; ++j){
C[i][j] = 0;
for (k = 0; k < N; ++k){
C[i][j] += A[i][k] * B[k][j];
}
}
}
但是这里有一个局限,如果假设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]。
故进一步的并行算法为:
#prama omp parallel for num_threads(NUM)
for (z = 0; z < M * P; ++z){
i = z / P;
j = z % P;
C[i][j] = 0;
for (k = 0; k < N; ++k){
C[i][j] += A[i][k] * B[k][j];
}
}
性能优化。
看最里面一层的计算
for (k = 0; k < N; ++k){
C[i][j] += A[i][k] * B[k][j];
由于内存中二维数组是以行优先进行存储的,因此B[k][j]存在严重的cache命中率问题,解决这个问题的方法是也将B进行一次沿对角线进行翻转,使得最里面的计算变成
for (k = 0; k < N; ++k){
C[i][j] += A[i][k] * B[j][k];
另外一点需要注意的就是C[i][j] += A[i][k] * B[j][k];计算时的伪共享问题。