我想知道是否有人可以向我展示如何有效地使用循环平铺/循环阻塞来进行大型密集矩阵乘法。我正在用 1000x1000 矩阵做C = AB 。我已经按照 Wikipedia 上的示例进行循环平铺，但使用平铺得到的结果比不使用时更差。
http://en.wikipedia.org/wiki/Loop_tiling
我在下面提供了一些代码。由于缓存未命中，天真的方法非常慢。transpose 方法在缓冲区中创建B的转置。此方法给出最快的结果（矩阵乘法为 O(n^3)，转置为 O(n^2)，因此转置速度至少快 1000 倍）。没有阻塞的wiki方法也很快，不需要缓冲区。阻塞方法较慢。阻塞的另一个问题是它必须多次更新块。这对线程/OpenMP 来说是一个挑战，因为如果不小心，它可能会导致竞争条件。
我应该指出，使用 AVX 和转置方法的修改，我得到的结果比 Eigen 更快。但是，我使用 SSE 的结果比 Eigen 慢一些，所以我认为我可以更好地使用缓存。
编辑：我想我知道我想做什么。它来自 1969 年的 Cannon 算法。
http://en.wikipedia.org/wiki/Matrix_multiplication#Communication-avoiding_and_distributed_algorithms
我需要进行块矩阵乘法并让每个线程处理C而不是A和B的子矩阵。例如，如果我将矩阵分成四个块。我会做：
+-+ +-+ +-+ +-+ +-+ +-+
| | | | | |
| C1 C2 | | A1 A2 | | B1 B2 |
| | = | | x | |
| C3 C4 | | A3 A4 | | B3 B4 |
| | | | | |
+-+ +-+ +-+ +-+ +-+ +-+
thread 1: C1 = A1B1 + A2B3
thread 2: C2 = A1B2 + A2B4
thread 3: C3 = A3B1 + A4B3
thread 4: C4 = A3B2 + A4B4
这消除了竞争条件。我得考虑一下。
void matrix_mult_naive(const float*A , const float* B, float* C, const int N, const int M, const int K) {
#pragma omp parallel for
for(int i=0; i<N; i++) {
for(int j=0; j<K; j++) {
float tmp = 0;
for(int l=0; l<M; l++) {
tmp += A[M*i+l]*B[K*l+j];
}
C[K*i + j] = tmp;
}
}
}
void matrix_mult_transpose(const float*A , const float* B, float* C, const int N, const int M, const int K) {
float *B2 = (float*)aligned_malloc(M*K*sizeof(float), 32);
transpose(B, B2, M, K, 1);
#pragma omp parallel for
for(int i=0; i<N; i++) {
for(int j=0; j<K; j++) {
float tmp = 0;
for(int l=0; l<M; l++) {
tmp += A[M*i+l]*B2[M*j+l];
}
C[K*i + j] = tmp;
}
}
aligned_free(B2);
}
void matrix_mult_wiki(const float*A , const float* B, float* C, const int N, const int M, const int K) {
for(int i=0; i<N; i++) {
for(int j=0; j<K; j++) {
C[K*i + j] = 0;
}
}
#pragma omp parallel for
for(int i=0; i<N; i++) {
for(int l=0; l<M; l++) {
float a = A[M*i+l];
for(int j=0; j<K; j++) {
C[K*i + j] += a*B[K*l+j];
}
}
}
}
void matrix_mult_wiki_block(const float*A , const float* B, float* C, const int N, const int M, const int K) {
const int block_size = 8; //I have tried several different block sizes
for(int i=0; i<N; i++) {
for(int j=0; j<K; j++) {
C[K*i + j] = 0;
}
}
for(int l2=0; l2<M; l2+=block_size) {
for(int j2=0; j2<K; j2+=block_size) {
#pragma omp parallel for
for(int i=0; i<N; i++) {
for(int l=l2; l<min(M, l2+block_size); l++) {
for(int j=j2; j<min(K, j2+block_size); j++) {
C[K*i + j] += A[M*i+l]*B[K*l+j];
}
}
}
}
}
}