如何优化矩阵乘法 (matmul) 代码以在单个处理器内核上快速运行 [英] how to optimize matrix multiplication (matmul) code to run fast on a single processor core
问题描述
我正在研究并行编程概念并尝试在单核上优化矩阵乘法示例.到目前为止，我提出的最快的实现如下:
/* 此例程执行 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.1204CPU 上最先进的矩阵乘法实现使用
上图(原文来自这篇论文，直接在本教程) 说明了在 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 上优化矩阵乘法的更多信息，以下教程可能会有所帮助:
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).
The above figure (originally from this paper, directly used in this tutorial) illustrates the GotoBLAS algorithm as implemented in BLIS. Cache blocking parameters {MC, NC, KC} determine the submatrix sizes of Bp (KC × NC) and Ai (MC × KC), such that they fit in various caches. During the computation, row panels Bp are contiguously packed into buffer Bp to fit in the L3 cache. Blocks Ai are similarly packed into buffer Ai to fit in the L2 cache. Register block sizes {MR, NR} relate to submatrices in registers that contribute to C. In the micro-kernel (the inner most loop), a small MR × NR micro-tile of C is updated by pair of MR × KC and KC × NR slivers of Ai and Bp.
For the Strassen's algorithm with O(N^2.87) complexity, you might be interested in reading this paper. Other fast matrix multiplication algorithms with asymptotic complexity less than O(N^3) can be easily extended in this paper. There is a recent thesis about the practical fast matrix multiplication algorithms.
The following tutorials might be helpful if you want to learn more about how to optimize matrix multiplication on CPUs:
GEMM: From Pure C to SSE Optimized Micro Kernels
BLISlab: A sandbox for optimizing GEMM for CPU and ARM
A most updated document about how to optimize GEMM on CPUs (with AVX2/FMA) step by step can be downloaded here: https://github.com/ULAFF/LAFF-On-HPC/blob/master/LAFF-On-PfHP.pdf
A Massive Open Online Course to be offered on edX starting in June 2019 (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
这篇关于如何优化矩阵乘法 (matmul) 代码以在单个处理器内核上快速运行的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
