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

Question

我想知道是否有人可以告诉我如何有效地使用循环平铺/循环阻塞来进行大密集矩阵乘法.我正在使用1000x1000矩阵进行C = AB.我已经按照维基百科上的示例进行循环平铺,但是使用平铺比使用平铺更糟糕.

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

我在下面提供了一些代码.由于缓存未命中,天真的方法非常慢.转置方法在缓冲区中创建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 ^而非阿和乙.例如,如果我将矩阵分成四个块.我会做:

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

Answer 1

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

如果您必须选择 AVX 或阻止之一，那么使用 AVX 扩展（vfmadd###ps和haddps）仍然要快得多。64 / sizeof(float)考虑到您已经在测试块大小是否是== 16 个浮点数 == 两个 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 算法是更好的选择。

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