从硬件的视角看GEMM

从硬件的视角看GEMM
1. 分块乘法的内存层次架构
分块矩阵乘法,如图6-28所示, 通过将矩阵分块拆分,能够在处理器的Cache和寄存器内存放进行快速计算,计算完成后写回主存。
图6-28 分块矩阵乘法
首先,所有的数据都在主内存中,如图6-29所示。
图6-29 所有的数据都在主内存中
然后,在分块加载到L2Cache中, 完成分块子任务,如图6-30所示。
图6-30 在分块加载到L2Cache中, 完成分块子任务
在进入计算核内部时,将向量Block进一步拆分成更小的Tile进行计算,如图6-32所示。
图6-32 在进入计算核内部时,将向量Block进一步拆分成更小的Tile进行计算
整个计算过程的访存层次化结构,如图6-33所示。
图6-33 整个计算过程的访存层次化结构
这样就对英伟达的矩阵乘法流程清楚了,通过多次将矩阵逐渐分为更小的Tile进行计算,如图6-34所示。
 
图6-34 通过多次将矩阵逐渐分为更小的Tile进行计算
2. 内存布局
在计算的过程中, 为了保证访存的连续性,A矩阵需要按行排序, 而B矩阵需要按列排序,如图6-35所示。
图6-35 为了保证访存的连续性,A矩阵需要按行排序, 而B矩阵需要按列排序
在分块的内部进行块乘法时,访存顺序变为列优先/行优先的方式, 因此矩阵的布局变成了一种Z字排列, 如图6-36所示。
 
图6-36 矩阵的布局变成了一种Z字排列
6.3.4 通用矩阵乘法执行
使用两个手工实现的纯粹GEMM和分块GEMM的例子来解释矩阵分块乘法的原理和性能影响, 可以看到性能差距接近53倍. 按照测试的A10 GPU峰值FP32算力31TFFLOPS来算, 最朴素的算法由于访存效率的问题, 浮点算力仅为峰值的1%。
# ./naive 
AveragePerformance  0.2336 Tflops

# ./block 
AveragePerformance  10.7669 Tflops
1. 纯粹GEMM
最简单的矩阵乘法如下:
#define OFFSET(row, col, stride) ((row) * (stride) + (col))

__global__ void basic_gemm(
    float *  A, float *  B, float *  C,
    const int M, const int N, const int K) {

    int _x = blockIdx.x * blockDim.x + threadIdx.x;
    int _y = blockIdx.y * blockDim.y + threadIdx.y;
    if (_x < M && _y < N) {
        float sum = 0.0;
        for (int k = 0; k < K; k++) {
            sum +=A[OFFSET(_x, k, K)] *  B[OFFSET(k , _y, N)];
        }
        C[OFFSET(_x, _y, N)] = sum;
    }
}
在A10上测试其FLOS大概仅有233GFlops。
int main() {
    const int M = 4096;
    const int K = 1024;
    const int N = 4096;
    const int ITER = 100;

    dim3 gridDim(ceil(M/32), ceil(N/32), 1);
    dim3 blockDim(32, 32, 1);

    float *d_a, *d_b, *d_c ;
    cudaMalloc(&d_a, M * K * sizeof(float));
    cudaMalloc(&d_b, K * N * sizeof(float));
    cudaMalloc(&d_c, M * N * sizeof(float));

    cudaEvent_t start, end;
    cudaEventCreate(&start);
    cudaEventCreate(&end);

    cudaEventRecord(start);
    for (int i = 0; i < ITER; i++)
        basic_gemm<<<gridDim, blockDim>>>(d_a, d_b, d_c, M, N, K);
    cudaEventRecord(end);
    cudaEventSynchronize(end);

    float msec;
    cudaEventElapsedTime(&msec, start, end);
    
    long workload =  long(M) * N * K * 2 * ITER;
    double avg_Tflops = ((double)workload / 1e12 ) / (double(msec)/ 1e3);
    printf("AveragePerformance  %6.4lf Tflops\n",avg_Tflops);

    cudaFree(d_a);
    cudaFree(d_b);
    cudaFree(d_c);
}
2. 分块GEMM
相关的代码画了一个容易理解的示意图,如图6-37所示。
 
图6-37 分块通用矩阵乘法示例代码示意图
代码如下:
__global__ void block2d_gemm(const float *A, const float *B, float *C,
                            int M, int N, int K) {

  const int BM = 128;
  const int BN = 128;
  const int BK = 8;
  const int TM = 8;
  const int TN = 8;

  const uint cRow = blockIdx.y;
  const uint cCol = blockIdx.x;

  const uint totalResultsBlocktile = BM * BN;
//线程负责计算块中的TM*TN元素
    const uint numThreadsBlocktile = totalResultsBlocktile / (TM * TN);

//BN/TN是跨越一列的线程数
  const int threadCol = threadIdx.x % (BN / TN);
  const int threadRow = threadIdx.x / (BN / TN);
 //在smem中为当前块分配空间
   __shared__ float As[BM * BK];
  __shared__ float Bs[BK * BN];
 //将块图块移动到A行和B列的开头
   A += cRow * BM * K;
  B += cCol * BN;
  C += cRow * BM * N + cCol * BN;
  //计算此线程将加载到SMEM中的索引
   const uint innerRowA = threadIdx.x / BK;
  const uint innerColA = threadIdx.x % BK;
   //计算单个块在单个步骤中加载的As行数
    const uint strideA = numThreadsBlocktile / BK;
  const uint innerRowB = threadIdx.x / BN;
  const uint innerColB = threadIdx.x % BN;
  //As和Bs,希望每个负载都跨越整个列宽,以便更好地进行GMEM合并
//(与跨越整个行宽和跨列迭代相反)
const uint strideB = numThreadsBlocktile / BN;
//为注册表文件中的结果分配线程本地缓存
  float threadResults[TM * TN] = {0.0};
//为As和Bs注册缓存
  float regM[TM] = {0.0};
  float regN[TN] = {0.0};
  //最外层的分块平铺环
  for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
//填充SMEM缓存
    for (uint loadOffset = 0; loadOffset < BM; loadOffset += strideA) {
      As[(innerRowA + loadOffset) * BK + innerColA] =
          A[(innerRowA + loadOffset) * K + innerColA];
    }
    for (uint loadOffset = 0; loadOffset < BK; loadOffset += strideB) {
      Bs[(innerRowB + loadOffset) * BN + innerColB] =
          B[(innerRowB + loadOffset) * N + innerColB];
    }
    __syncthreads();
  //超前分块
    A += BK;     //将BK列向右移动
    B += BK * N; //向下移动BK行

    // 计算每个线程的结果
    for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
// 进入内存的块
      for (uint i = 0; i < TM; ++i) {
        regM[i] = As[(threadRow * TM + i) * BK + dotIdx];
      }
      for (uint i = 0; i < TN; ++i) {
        regN[i] = Bs[dotIdx * BN + threadCol * TN + i];
      }
      for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
        for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
          threadResults[resIdxM * TN + resIdxN] +=
              regM[resIdxM] * regN[resIdxN];
        }
      }
    }
    __syncthreads();
  }

  // 写出结果
  for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
    for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
      C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN] =
         threadResults[resIdxM * TN + resIdxN] +
         C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN];
    }
  }
}
 
6.4 张量 矩阵乘法优化
在SIMT架构下, 不使用TensorCore进行矩阵乘法,计算所需要的访存相关的优化。通过逐步迭代优化,深入理解GPU的性能特征和内存访问优化。
测试环境为一块A10 GPU,  驱动版本: 550.54.15, CUDA版本: 12.4 . 矩阵M=N=K=4092,见表6-5。
表6-5 cuBLAS调用,在每种大小下调用的内核

内核

GFLOPs/s

相对于cuBLAS的性能

0:cuBLAS

14765.9

100.0%

1:Naive

588.5

3.9%

2:GMEM Coalescing

1165.7

7.9%

3:SMEM Caching

2166.7

14.6%

4:1D Blocktiling

6082.0

41.2%

5:2D Blocktiling

11279.0

76.4%

6:Vectorized Mem Access

12861.4

87.1%

7:WarpTiling

14766.3

100.0%

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.hqwc.cn/news/809367.html

如若内容造成侵权/违法违规/事实不符,请联系编程知识网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

电表行业强力推荐KY8023S6对标平替BL8023双向继电器驱动芯片

双向继电器驱动芯片KY8023S6对标平替BL8023。 产品特点 SOT23-6L和SOP-8 两种封装 静态功耗电流<1uA 输入高低转换电平为1.3V(Is=1A)左右,兼容多种单片机 工作电压5-40V,推荐工作电压5-30V 输入INA和INB端自带115k下拉电阻,一般可省略外围下拉电阻 集成高速续流二极管,…

重启后,idea 依赖爆红

没改动代码,但无法运行java项目 java: -source 1.5 中不支持……表达式/运算符java: 错误: 不支持发行版本 5解决:https://www.cnblogs.com/liu-han/p/16105953.html 文件→项目结构,SDK和项目/模块 语言级别的版本要一致SDK: 1.8 项目语言级别:SDK默认值 模块语言级别:SD…

[亲测]君临九州+单机安装教程+无需虚拟机+GM后台

今天给大家带来一款单机游戏的架设:君临九州,版本号:1.0.12.30。 另外:本人承接各种游戏架设(单机+联网) 本人为了学习和研究软件内含的设计思想和原理,带了架设教程仅供娱乐。 教程是本人亲自搭建成功的,绝对是完整可运行的,踩过的坑都给你们填上了。如果你是小白也没…

智慧园区管理系统原型

智慧园区管理系统的构建是一个复杂而系统的工程,它融合了信息化、AI、物联网等多种先进技术,旨在提升园区的管理效率、服务质量以及企业运营效率。 一、明确系统目标和需求 需求收集与分析:首先,需要对园区的实际需求进行全面分析,包括园区类型(如产业园区、办公园区、住…

POI导出excel文件加水印

百分百能用,我用的POI版本是5.2.3,效果如下import lombok.extern.slf4j.Slf4j; import org.apache.poi.openxml4j.opc.PackagePartName; import org.apache.poi.openxml4j.opc.PackageRelationship; import org.apache.poi.openxml4j.opc.TargetMode; import org.apache.poi.xs…

FredNormer: 非平稳时间序列预测的频域正则化方法

时间序列预测是一个具有挑战性的任务,尤其是在处理非平稳数据时。现有的基于正则化的方法虽然在解决分布偏移问题上取得了一定成功但仍存在局限性。这些方法主要在时间域进行操作,可能无法充分捕捉在频域中更明显的动态模式,从而导致次优的结果。 FredNormer论文的研究目的主要…

江苏省第二届数据安全技术应用职业技能竞赛初赛WP

一、数据安全解题赛1、ds_0602解题思路题目让我们获取加密文件中的原始数据,解密后提交第六行第二列数据,下载附件,发现里面有两个文件,其中一个是“.enc”结尾,那这里我们得先简单了解一下“.enc”结尾的是什么类型的文件。简单来说“.enc”结尾的文件通常是经过加密的文…

java之使用CompletableFuture入门2

Java 17 -序章 本文介绍用过的 allOf、anyOf 函数的用法。allOf 函数原型两点: 1、没有返回值。 2、参数 cfs 中 任何一个 都不能是 null。anyOf 函数原型两点: 1、有返回值,为 Object。 2、参数 cfs 中 任何一个 都不能是 null。allOf 测试意图: 多个任务正常执行。ben发布…

VMware Aria Operations for Logs 8.18 发布,新增功能概览

VMware Aria Operations for Logs 8.18 发布,新增功能概览VMware Aria Operations for Logs 8.18 - 集中式日志管理 请访问原文链接:https://sysin.org/blog/vmware-aria-operations-for-logs/,查看最新版。原创作品,转载请保留出处。 作者主页:sysin.org集中式日志管理 V…

VMware Aria Operations for Networks 6.13 发布,新增功能概览

VMware Aria Operations for Networks 6.13 发布,新增功能概览VMware Aria Operations for Networks 6.13 - 网络和应用监控工具 请访问原文链接:https://sysin.org/blog/vmware-aria-operations-for-networks/,查看最新版。原创作品,转载请保留出处。 作者主页:sysin.org…

读数据工程之道:设计和构建健壮的数据系统01数据工程概述

数据工程概述1. 数据工程 1.1. 自从公司开始使用数据做事,数据工程就以某种形式存在了1.1.1. 预测性分析、描述性分析和报告1.2. 数据工程师获取数据、存储数据,并准备数据供数据科学家、分析师和其他人使用 1.3. 数据工程是系统和流程的开发、实施和维护,这些系统和流程接收…

安装socks5的一次尝试

1. 下载并自动配置socks5sudo wget https://ap-guangzhou-1257892306.cos.ap-guangzhou.myqcloud.com/asi/httpsocks5.sh && sh httpsocks5.sh 执行下载脚本 wget —no-check-certificate https://raw.github.com/Lozy/danted/master/install.sh -O install.sh执行安装…