caching - block 矩阵矩阵乘法的最佳 block 大小值

Question

我想使用以下 C 代码进行块矩阵-矩阵乘法。在这种方法中，将大小为 BLOCK_SIZE 的块加载到最快的缓存中，以减少计算过程中的内存流量。

void bMMikj(double **A , double **B , double ** C , int m, int n , int p , int BLOCK_SIZE){

   int i, j , jj, k , kk ;
   register double jjTempMin = 0.0 , kkTempMin = 0.0;

   for (jj=0; jj<n; jj+= BLOCK_SIZE) {
       jjTempMin = min(jj+ BLOCK_SIZE,n); 
       for (kk=0; kk<n; kk+= BLOCK_SIZE) {
           kkTempMin = min(kk+ BLOCK_SIZE,n); 
           for (i=0; i<n; i++) {
               for (k = kk ; k < kkTempMin ; k++) {
                   for (j=jj; j < jjTempMin; j++) {
                      C[i][j]  +=  A[i][k] * B[k][j];
                   }
               }
           }
      }
   }
}

我搜索了 BLOCK_SIZE 的最佳合适值我发现 BLOCK_SIZE <= sqrt( M_fast / 3 )和 M_fast这是 L1 缓存。

在我的计算机中，我有两个 L1 缓存，如图 here与 lstopo 工具。

下面，我使用启发式方法，例如从 BLOCK_SIZE 开始的 4并将值增加 8为 1000次，具有不同的矩阵大小值。

跳跃获得最佳MFLOPS(或最少乘法时间)值和相应的 BLOCK_SIZE value 将是最合适的值。

这是用于测试的代码:

int BLOCK_SIZE = 4;
int m , n , p;
m = n = p = 1024; /* This value is also changed
                     and all the matrices are square, for simplicity
                     */
for(int i=0;i< 1000; i++ , BLOCK_SIZE += 8) {
    # aClock.start();
    test_bMMikj(A , B ,  C , loc_n , loc_n , loc_n ,BLOCK_SIZE);
    # aClock.stop();
}

测试给了我每个矩阵大小的不同值并且不同意公式。计算机型号是'Intel® Core™ i5-3320M CPU @ 2.60GHz × 4' with 3.8GiB，这里是 Intel specification

另一个问题:
如果我有两个 L1 缓存，就像我在 here 中所拥有的一样，我应该考虑 BLOCK_SIZE关于其中之一还是两者的总和？

Answer 1

1.块矩阵乘法:
这个想法是通过重用当前存储在缓存中的数据块来最大限度地利用时间和空间局部性。您的代码不正确，因为它只包含 5 个循环；对于块应该有 6 个，例如:

for(int ii=0; ii<N; ii+=stride)
{
    for(int jj=0; jj<N; jj+=stride)
    {
        for(int kk=0; kk<N; kk+=stride)
        {
            for(int i=ii; i<ii+stride; ++i)
            {
                for(int j=jj; j<jj+stride; ++j)
                {
                    for(int k=kk; k<kk+stride; ++k) C[i][j] += A[i][k]*B[k][j];
                }
            }               
        }
    }
}

为简单起见，最初将 N 和 stride 都保留为 2 的幂。
ijk 模式不是最优的，你应该选择 kij 或 ikj，详细信息 here .不同的访问模式有不同的表现，你应该尝试ijk的所有排列。

2.块/步幅大小:
一般说来，最快的缓存 (L1) 应该能够容纳 3 个数据块(步幅*步幅)，以便在矩阵乘法的情况下获得最佳性能，但自己尝试并找到它总是好的。
将 stride 增加 8 可能不是一个好主意，请尝试将其保持为 2 的增加幂，因为大多数块大小都是以这种方式调整大小的。
而且您应该只查看数据缓存(L1d)，在您的情况下为 32KB。