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

本文介绍了如何优化矩阵乘法 (matmul) 代码以在单个处理器内核上快速运行的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着跟版网的小编来一起学习吧！

问题描述

限时送ChatGPT账号..

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

/* 此例程执行 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).





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:

How to Optimize GEMM Wiki

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) 代码以在单个处理器内核上快速运行的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持跟版网！
                本站部分内容来源互联网,如果有图片或者内容侵犯了您的权益，请联系我们，我们会在确认后第一时间进行删除！