如何优化矩阵乘法 (matmul) 代码以在单个处理器内

问题描述

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

/* 此例程执行 dgemm 操作* C := C + A * B* 其中 A、B 和 C 是以列优先格式存储的 lda-by-lda 矩阵.* 退出时，A 和 B 保持它们的输入值.*/void square_dgemm (int n, double* A, double* B, double* C){/* 对于 A 的每一行 i */for (int i = 0; i 

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

login4.stampede(72)$ tail -f job-naive.stdout大小:480 Mflop/s:1818.89 百分比:18.95大小:511 Mflop/s:2291.73 百分比:23.87大小:512 Mflop/s:937.061 百分比:9.76大小:639 Mflop/s:293.434 百分比:3.06大小:640 Mflop/s:270.238 百分比:2.81大小:767 Mflop/s:240.209 百分比:2.50大小:768 Mflop/s:242.118 百分比:2.52大小:769 Mflop/s:240.173 百分比:2.50峰值的平均百分比 = 22.0802等级 = 33.1204

CPU 上最先进的矩阵乘法实现使用

上图(原文来自这篇论文，直接在本教程) 说明了在 BLIS.缓存阻塞参数{MC, NC, KC}确定Bp (KC × NC) 和 Ai (MC × KC) 的子矩阵大小，以便它们适合各种缓存.在计算过程中，行面板 Bp连续打包到缓冲区 Bp 中以适合 L3 缓存.块 Ai 类似地打包到缓冲区 Ai以适合 L2 缓存.寄存器块大小 {MR, NR} 与寄存器中对 C 有贡献的子矩阵有关.在微内核(最内部的循环)中，C 的一个小的 MR × NR 微瓦片由一对 MR × KC 和 KC 更新× Ai 和 Bp 的 NR 条带.

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

如果您想了解有关如何在 CPU 上优化矩阵乘法的更多信息，以下教程可能会有所帮助:

如何优化 GEMM Wiki

GEMM:从纯 C 到 SSE 优化的微内核

BLISlab:针对 CPU 和 ARM 优化 GEMM 的沙箱

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

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

I am working on parallel programming concepts and trying to optimize matrix multiplication example on single core. The fastest implementation I came up so far is the following:

/* 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;
    }
}

The results are like below. how to reduce the loops and increase the performance

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

The state-of-the-art implementation of matrix multiplication on CPUs uses GotoBLAS algorithm. Basically the loops are organized in the following order:

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

A key insight underlying modern high-performance implementations of matrix multiplication is to organize the computations by partitioning the operands into blocks for temporal locality (3 outer most loops), and to pack (copy) such blocks into contiguous buffers that fit into various levels of memory for spatial locality (3 inner most loops).