块矩阵乘法的最佳块大小值

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缓存，如此处所示，使用lstopo工具可以查看。 以下是使用启发式方法的示例，从BLOCK_SIZE为4开始，每次增加8，重复1000次，使用不同的矩阵大小。

希望得到最佳MFLOPS（或最少的乘法时间）值和相应的BLOCK_SIZE值，这将是最适合的值。

下面是测试代码：

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”，内存为3.8GiB，此处是Intel规格说明。

另外一个问题：
如果我有两个L1缓存，像我在这里看到的那样，我应该将BLOCK_SIZE考虑为其中一个还是两者之和？

- Nawal

---

使用 lstopo 检查 CPU 架构（NUMA 缓存）并加 1。除非您披露编译详细信息，否则没有人能告诉您更多。因此，在物理核心 P＃0 上为代码执行定位的线程将无法从属于 P＃1 物理核心的 L1d 中的任何数据中受益，因此关于“共享”存储的猜测仅从 L3 缓存开始有效（实际上不超过约 3 MB）。还要始终检查实际的 CPU 缓存关联性、缓存行大小和所有细节，以确定使用缓存预取掩盖 DRAM 访问延迟的机会。 - user3666197

谢谢。你的意思是我只需要考虑L1d和L2d吗？你能告诉我如何检查实际的CPU缓存关联性和缓存行大小吗？ - Nawal

在x86上，CPUID可以在运行时提供有关此信息。 https://en.wikipedia.org/wiki/CPUID#EAX.3D80000005h:_L1_Cache_and_TLB_Identifiers. （从您的缓存大小/层次结构来看，我认为您正在使用双核（具有超线程）英特尔CPU，例如Haswell i3台式机或i3 / i5笔记本电脑。）如果您希望其运行速度快，并且进行缓存阻塞，您需要让编译器自动矢量化它（或使用SSE2，AVX和理想情况下的FMA手动矢量化）。通过良好的缓存阻塞，matmul可以更多或更少地饱和乘法/加法或FMA吞吐量，即使使用32字节向量也是如此。 - Peter Cordes

是的，谢谢。我添加了计算机型号细节。 - Nawal

不确定这是否有帮助，但我听说寄存器双精度没有用处，这是真的吗？ - Aditya P

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和步幅都设置为2的次方。 ijk模式不是最优的，您应该选择kij或ikj，这里有详细信息 here。不同的访问模式具有不同的性能，您应该尝试所有的ijk排列。

2. 块/步长大小： 通常认为，您最快的缓存（L1）应该能够容纳3个块（stride*stride）的数据以获得矩阵乘法的最佳性能，但是进行实验并自行找到它总是很好的。 将步幅增加8可能不是一个好主意，尝试将其保持为2的次方的增量，因为大多数块大小都是以此方式确定的。 对于您的情况，您只需要查看数据缓存（L1d），它的大小为32KB。

- Kshitij Shukla