c++ - 如何优化矩阵乘法 (matmul) 代码以在单个处理器内核上快速运行

Question

我正在研究并行编程概念并尝试优化单核上的矩阵乘法示例。到目前为止，我想出的最快的实现如下:

/* This routine performs a dgemm operation
 *  C := C + A * B
 * where A, B, and C are lda-by-lda matrices stored in column-major format.
 * On exit, A and B maintain their input values. */    
void square_dgemm (int n, double* A, double* B, double* C)
{
  /* For each row i of A */
  for (int i = 0; i < n; ++i)
    /* For each column j of B */
    for (int j = 0; j < n; ++j) 
    {
      /* Compute C(i,j) */
      double cij = C[i+j*n];
      for( int k = 0; k < n; k++ )
          cij += A[i+k*n] * B[k+j*n];
      C[i+j*n] = cij;
    }
}

结果如下。如何减少循环并提高性能

login4.stampede(72)$ tail -f job-naive.stdout
Size: 480       Mflop/s:  1818.89       Percentage: 18.95
Size: 511       Mflop/s:  2291.73       Percentage: 23.87
Size: 512       Mflop/s:  937.061       Percentage:  9.76
Size: 639       Mflop/s:  293.434       Percentage:  3.06
Size: 640       Mflop/s:  270.238       Percentage:  2.81
Size: 767       Mflop/s:  240.209       Percentage:  2.50
Size: 768       Mflop/s:  242.118       Percentage:  2.52
Size: 769       Mflop/s:  240.173       Percentage:  2.50
Average percentage of Peak = 22.0802
Grade = 33.1204

Answer 1

CPU 上最先进的矩阵乘法实现使用 GotoBLAS算法。基本上，循环按以下顺序组织:

Loop5 for jc = 0 to N-1 in steps of NC
Loop4   for kc = 0 to K-1 in steps of KC
          //Pack KCxNC block of B
Loop3     for ic = 0 to M-1 in steps of MC
            //Pack MCxKC block of A
//--------------------Macro Kernel------------
Loop2       for jr = 0 to NC-1 in steps of NR
Loop1         for ir = 0 to MC-1 in steps of MR
//--------------------Micro Kernel------------
Loop0           for k = 0 to KC-1 in steps of 1
                //update MRxNR block of C matrix

矩阵乘法现代高性能实现的一个关键见解是通过将操作数划分为时间局部性的 block (3 个最外层循环)来组织计算，并打包(复制)这些 block 放入适合空间局部性的不同内存级别的连续缓冲区(3 个最内层循环)。

上图(来源于this paper，直接用于this tutorial)说明了BLIS中实现的GotoBLAS算法。 .缓存阻塞参数 {MC, NC, KC} 确定Bp (KC × NC) 和 Ai (MC × KC) 的子矩阵大小，以便它们适合各种缓存。在计算过程中，行面板 Bp被连续打包到缓冲区 Bp 中以适应 L3 缓存。 block Ai 被类似地打包到缓冲区 Ai 中以适应二级缓存。寄存器 block 大小 {MR, NR} 与寄存器中对 C 有贡献的子矩阵有关。在微内核(最内层循环)中，C 的一个小 MR × NR 微 block 由 MR × KC 和 KC 对更新× Ai 和 Bp 的 NR 条。

对于 O(N^2.87) 复杂度的 Strassen 算法，您可能有兴趣阅读 this paper .其他渐近复杂度小于 O(N^3) 的快速矩阵乘法算法可以在 this paper 中轻松扩展.有一个 recent thesis关于实用的快速矩阵乘法算法。

如果您想详细了解如何在 CPU 上优化矩阵乘法，以下教程可能会对您有所帮助:

How to Optimize GEMM Wiki

GEMM: From Pure C to SSE Optimized Micro Kernels

BLISlab: A sandbox for optimizing GEMM for CPU and ARM

有关如何在 CPU 上(使用 AVX2/FMA)逐步优化 GEMM 的最新文档可在此处下载: https://github.com/ULAFF/LAFF-On-HPC/blob/master/LAFF-On-PfHP.pdf

将于 2019 年 6 月开始在 edX 上提供大规模开放式在线类(class)(LAFF-On Programming for High Performance): https://github.com/ULAFF/LAFF-On-HPC http://www.cs.utexas.edu/users/flame/laff/pfhp/LAFF-On-PfHP.html