gpt4 book ai didi

c - 用于大型密集矩阵乘法的循环平铺/阻塞

转载 作者:太空狗 更新时间:2023-10-29 16:32:45 27 4
gpt4 key购买 nike

我想知道是否有人可以告诉我如何有效地使用循环平铺/循环阻塞进行大型密集矩阵乘法。我正在用 1000x1000 矩阵做C = AB。我按照 Wikipedia 上的循环平铺示例进行操作,但使用平铺得到的结果比不使用时更差。

http://en.wikipedia.org/wiki/Loop_tiling

http://software.intel.com/en-us/articles/how-to-use-loop-blocking-to-optimize-memory-use-on-32-bit-intel-architecture

我在下面提供了一些代码。由于缓存未命中,天真的方法非常慢。 transpose 方法在缓冲区中创建 B 的转置。这种方法给出最快的结果(矩阵乘法为 O(n^3),转置为 O(n^2),因此转置至少快 1000 倍)。没有阻塞的 wiki 方法也很快,不需要缓冲区。阻塞方法较慢。阻塞的另一个问题是它必须多次更新 block 。这对线程/OpenMP 来说是一个挑战,因为如果不小心,它可能会导致竞争条件。

我应该指出,使用 AVX 并修改转置方法可以比 Eigen 更快地获得结果。但是,我使用 SSE 的结果比 Eigen 慢一点,所以我认为我可以更好地使用缓存。

编辑:我想我知道我想做什么。它来自于 1969 年的 Cannon 算法。
http://en.wikipedia.org/wiki/Matrix_multiplication#Communication-avoiding_and_distributed_algorithms

我需要进行 block 矩阵乘法,并让每个线程处理 C 而不是 A 的子矩阵。 strong> 和 B。例如,如果我将矩阵分成四个 block 。我会这样做:

+-+      +-+     +-+      +-+   +-+      +-+
| | | | | |
| 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];
}
}
}
}
}
}

最佳答案

我得到的最佳结果是再添加一个 for 循环来阻塞 N,然后重新排列循环。我还提升了循环不变代码,但编译器的优化器应该会自动执行此操作。 block 大小应该是缓存行大小除以 sizeof(float)。这比转置方法快 50%。

如果您只能选择 AVX 或阻塞之一,使用 AVX 扩展(vfmadd###pshaddps)仍然会快得多。考虑到您已经在测试 block 大小是否是 64/sizeof(float) == 16 个 float == 两个 256 位 AVX 寄存器的倍数,因此最好和直接使用两者。

  • 转置:1,816,522 个刻度
  • 平铺:892,431(快 51%)
  • AVX+平铺:230,512(快 87%)

平铺:

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 = 64 / sizeof(float); // 64 = common cache line size
for(int i=0; i<N; i++) {
for(int j=0; j<K; j++) {
C[K*i + j] = 0;
}
}
for (int i0 = 0; i0 < N; i0 += block_size) {
int imax = i0 + block_size > N ? N : i0 + block_size;

for (int j0 = 0; j0 < M; j0 += block_size) {
int jmax = j0 + block_size > M ? M : j0 + block_size;

for (int k0 = 0; k0 < K; k0 += block_size) {
int kmax = k0 + block_size > K ? K : k0 + block_size;

for (int j1 = j0; j1 < jmax; ++j1) {
int sj = M * j1;

for (int i1 = i0; i1 < imax; ++i1) {
int mi = M * i1;
int ki = K * i1;
int kij = ki + j1;

for (int k1 = k0; k1 < kmax; ++k1) {
C[kij] += A[mi + k1] * B[sj + k1];
}
}
}
}
}
}
}

至于 Cannon 引用,SUMMA algorithm是一个更好的遵循。


万一其他人正在优化高瘦乘法({~1e9 x 50} x {50 x 50},我是如何在这里结束的),转置方法在性能上几乎与阻塞方法相同 n=18( float )。 n=18 是一种病态情况(比 17 或 19 更糟),我不太清楚导致这种情况的缓存访问模式。所有较大的 n 都通过阻塞方法得到改进。

关于c - 用于大型密集矩阵乘法的循环平铺/阻塞,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/15829223/

27 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com