CUDA编程之GEMM优化

置顶 已于 2022-11-04 20:22:04 修改
阅读量2.5k 收藏 26
点赞数 3
分类专栏： 高性能计算 深度学习 文章标签： 性能优化 深度学习
于 2022-10-17 16:18:32 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qianqing13579/article/details/127359866
版权

前言

最近由于工作需要，研究了一下CUDA编程中的GEMM的优化，主要是学习了GEMM优化的常用方法，同时自己也利用了常用优化方法动手实现了一遍GEMM优化。学习过程中主要参考了CUTLASS官方博客，网上也有中文翻译版本，里面有些地方翻译的可能不是很准确，在阅读中文版本的时候最好能对照原文看一下，在学习过程中还参考了网上的其他一些资料：

  1. CUDA SGEMM矩阵乘法优化笔记——从入门到cublas
  2. 深入浅出GPU优化系列：GEMM优化（一）

这篇文章主要想谈谈自己对GEMM优化的一些理解。

文章目录

Matrix类的设计

由于本文的所有示例代码都是基于Matrix类来设计的，所以这里有必要先讲一下Matrix类的设计。

为什么要设计Matrix

在自己手动实现GEMM优化的时候，发现一个问题：由于GEMM优化的时候需要对block和thread分块，这样就会涉及到大量子块以及子块元素索引的计算，很容易出错，为了简化索引的计算和提高可读性，设计了Matrix这个类型。

设计思想

在设计Matrix类的时候，我参考了Pytorch的Tensor的设计思想，Pytorch的作者在一篇博客中详细介绍了Tensor的设计思想，我们都知道Pytorch中的Tensor类型不保存实际的数据，也就意味着Tensor不负责实际数据的分配和释放，真实数据的分配和释放是由Storage这个类型来操作的，Tensor只是实际数据的一个逻辑视图(view)。下图来自那篇博客，图中很清楚的说明了Tensor和Storage的关系。
在这里插入图片描述
基于Pytorch的设计思想，Matrix类也被设计成了Tensor一样的机制，也就是Matrix不保存实际的内存数据，Matrix只是实际数据的一个视图(view)。Matrix除了拥有一个数据指针外，还具有两个重要的属性：size和stride。通过这两个属性，我们就可以很方便的访问矩阵数据了。下面我们看一下Matrix的视图机制是如何实现数据访问的。
下图表示内存中的一个4行6列的二维数组，该数组按照行主序的方式在内存中连续存储(与C语言中的数组一致)，所以在列这个维度上步长为1，在行这个维度上的步长为6。
在这里插入图片描述
现在我们要从该数组中获取一个4行3列的一个子数组，如下图所示：
在这里插入图片描述
下面我们可以利用视图完成该操作。
在这里插入图片描述
我们让视图的数据指针指向原始数据的第3个元素（原始数组下标为2的元素），然后设置视图的size为[4,3]，stride为[6,1]。如果我们要访问该子矩阵索引为(1,1)的元素（下图红色元素），则其内存索引为二维索引与stride的内积，即16+11=7，通过将data pointer移动7次就可以访问到该元素了。
在这里插入图片描述
通过上面对视图的介绍，我们可以看出通过视图机制，我们可以不用拷贝子数组的数据就可以很方便的访问该子数组。

代码实现

下面看一下Matrix的具体代码实现

template<typename T>
class Matrix
{
public:
    __device__ __host__ Matrix() = default;
    __device__ __host__ Matrix(const Matrix &) = default;
    __device__ __host__ Matrix& operator=(const Matrix &) = default;
    __device__ __host__ Matrix(T *_data,int _rows,int _cols,int _strideOfRow,int _strideOfCol):
                                    data(_data),
                                    rows(_rows),
                                    cols(_cols),
                                    strideOfRow(_strideOfRow),
                                    strideOfCol(_strideOfCol){}

    // 返回该矩阵所有字节数
    constexpr __device__ __host__ int GetNumberOfBytes() const
    {
        return rows*cols*sizeof(T);
    }

    // 返回该矩阵元素个数
    constexpr __device__ __host__ int GetNumberOfElements() const
    {
        return rows*cols;
    }

    // 访问某个元素，该元素的索引为二维逻辑索引：(rowIndex,colIndex)
    __device__ __host__ float &operator()(int rowIndex,int colIndex)
    {
        // 计算内存索引
        int memoryIndex=rowIndex*strideOfRow+colIndex*strideOfCol;

        return data[memoryIndex];
    }

    // 访问某个元素，该元素的索引为一维逻辑索引：(Index)
    __device__ __host__ float &operator()(int index)
    {
        // 转换为二维逻辑索引
        int colIndex=index%cols;
        int rowIndex=index/cols;

        // 计算内存索引
        int memoryIndex=rowIndex*strideOfRow+colIndex*strideOfCol;

        return data[memoryIndex];
    }



public:
    T *data = nullptr; // 数据指针
    int rows = 0;// 矩阵的行数
    int cols = 0;// 矩阵的列数
    int strideOfRow = 0;// 行步长
    int strideOfCol = 0;// 列步长

};
  1. 由于该类型需要同时在host端和device端使用，所以里面的成员都需要设置为__device__ host
  2. 由于该类型是一个视图类，所以拷贝构造和拷贝赋值直接使用默认的就可以了，不用自定义
  3. 这里需要区分一下两个概念：逻辑索引和内存索引，通常我们访问矩阵的数据的时候都是使用的逻辑索引，比如a[1][1],这里的索引(1,1)就是逻辑索引，而计算机在访问实际数据的时候，使用的是内存索引，逻辑索引(1,1)会被转换为实际的内存索引，比如上面的示例中，逻辑索引(1,1)被转换为了内存索引7

有了Matrix这个类，我们就可以方便的访问子矩阵的元素了。前面提到的示例可以通过下面的代码来实现，其中矩阵b就表示a中黄色部分的子矩阵。


Matrix<float> a;
Matrix<float> b(a.data+2,4,3,a.strideOfRow,a.strideOfCol );// a矩阵的一个4行3列的子矩阵
a(1,1)=0;// 访问索引为(1,1)的元素，实际访问的时候，会转换为内存索引7，然后通过data[7]访问到实际的元素

GEMM优化

介绍了Matrix的设计之后，下面我们来看一下GEMM优化。由于学习过程中主要参考了CUTLASS官方博客，所以本文的GEMM优化思路与这篇博客基本相同。

GEMM简介

在高性能计算领域，矩阵乘（GEMM）的优化是一个非常重要的课题。GEMM可以非常广泛地应用于航空航天、流体力
学及深度学习领域。GEMM的计算公式如下：
C = a l p h a A ∗ B + b e t a C C=alphaA*B+betaC C=alphaA∗B+betaC
在这里插入图片描述

其中A是一个M x K的矩阵，B是一个K x N的矩阵，C是一个M x N的矩阵，alpha和beta是标量。
为了简化计算，本文中设置alpha=1.0,beta=0。上述公式变为：
C = A ∗ B C=A*B C=A∗B

Baseline: NaiveGEMM

对于GEMM计算，我们很容易想到如下的实现：

__global__ void NaiveGEMM(Matrix<float> A,Matrix<float> B,Matrix<float> C)
{
    
    // 获取线程在网格内的索引
    int row = blockIdx.y * blockDim.y + threadIdx.y;// 行
    int col = blockIdx.x * blockDim.x + threadIdx.x;// 列

    // 每个线程计算矩阵C的一个元素
    if(row<C.rows&&col<C.cols)
    {
        float c = 0;
        for (int i = 0; i < A.cols; ++i)
        {
            c += A(row,i)*B(i,col);// 使用A的第row行乘以B的第col列
        }
        C(row,col) = c;
    }  
}

在这里插入图片描述

这种实现中，每个线程使用A的一行和B的一列计算矩阵C中的一个元素。这种实现存在如下问题：访存开销太大，A和B矩阵的全局内存被访问了多次，比如计算下图所示的C矩阵的一个绿色的3x3的子块，A的每一行被访问了3次(比如A的第一行在计算1,2,3子块的时候都被访问了一次)，B的每一列也被访问了3次，所以A矩阵被访问了N次，B矩阵被访问了M次。下面我们以这个实现为baseline来进行GEMM的优化。
在这里插入图片描述

减少访存策略1：Block分块

NaiveGEMM的主要问题在于访存开销太大，所以下面我们主要的优化点就是要减少访存，其实很容易想到的思路就是利用共享内存(Shared Memory)，基本思路就是将每个block对应的A矩阵和B矩阵的全局内存数据先读取到共享内存中，然后再计算，这样每次只需要读取一次全局内存就可以了，但是这里有个问题就是由于共享内存的容量是非常有限的，如果矩阵A和矩阵B规模比较大，则无法一次性全部加载到共享内存，所以我们需要分批加载。如下图所示，在计算C矩阵的绿色子块的时候，加载A和B矩阵的时候，沿着K维度每次只加载下图中1和2大小的子块。
在这里插入图片描述
所以GEMM优化的第一步就是Block分块：每个Block负责计算C矩阵的一个子块，每次计算子块的时候，沿着K维度分批将A和B加载到共享内存中计算乘累加。

与NaiveGEMM一样，这种实现方式每个线程还是负责block中一个元素的计算。通过Block分块，可以显著减少对全局内存的访问次数。假设每个Block的行大小为BM,列大小为BN,则优化后A矩阵被访问N/BN次，B矩阵被访问了M/BM次。

代码实现如下：

__global__ void BlockGEMM_V1(Matrix<float> A,Matrix<float> B,Matrix<float> C)
{
    // 注意命名不要与前面的宏定义重名
    const int BLOCK_M=16;// block的行数
    const int BLOCK_N=16;// block的列数
    const int BLOCK_K=16;

    // 沿着K维度循环加载一个block中对应的A和B的数据到共享内存
    float c=0.0;
    for(int i=0;i<A.cols/BLOCK_K;++i)
    {
        // 每个block对应的全局内存中的A,B子块，即创建全局内存中A,B的view
        Matrix<float> ASub(A.data+blockIdx.y*BLOCK_M*A.strideOfRow+i*BLOCK_K,BLOCK_M,BLOCK_K,A.strideOfRow,A.strideOfCol);
        Matrix<float> BSub(B.data+i*BLOCK_K*B.strideOfRow+blockIdx.x*BLOCK_N,BLOCK_K,BLOCK_N,B.strideOfRow,B.strideOfCol);

        // 将Asub,BSub加载到共享内存
        // 注意：这里需要将一维逻辑索引转换为多维逻辑索引：stardIndex->(stardIndex/cols,stardIndex%cols)
        __shared__ float A_Shared[BLOCK_M][BLOCK_K];
        __shared__ float B_Shared[BLOCK_K][BLOCK_N];
        int numberOfElementsPerThread=(BLOCK_K*BLOCK_M)/(blockDim.x*blockDim.y);// 每个线程需要读取多少数据
        int stardIndex=numberOfElementsPerThread*(threadIdx.y*blockDim.x+threadIdx.x);// stardIndex为每个线程读取的起始索引
        for(int threadIndex=0;threadIndex<numberOfElementsPerThread;++threadIndex)
        {
            int logicalIndex=stardIndex+threadIndex;
            A_Shared[logicalIndex/BLOCK_K][logicalIndex%BLOCK_K]=ASub(logicalIndex/BLOCK_K,logicalIndex%BLOCK_K);
            B_Shared[logicalIndex/BLOCK_N][logicalIndex%BLOCK_N]=BSub(logicalIndex/BLOCK_N,logicalIndex%BLOCK_N);
        }
        __syncthreads();

        // 每个thread计算A的一行和B的一列
        for(int k=0;k<BLOCK_K;++k)
        {
            c+=A_Shared[threadIdx.y][k]*B_Shared[k][threadIdx.x];
        }
        __syncthreads();

    }

    // 将每个线程计算好的结果写回到C矩阵
    // CSub为每个线程对应的全局内存的C矩阵子块，创建C矩阵的view
    Matrix<float> CSub(C.data+(blockIdx.y*BLOCK_M*C.strideOfRow+blockIdx.x*BLOCK_N),BLOCK_M,BLOCK_N,C.strideOfRow,C.strideOfCol);
    CSub(threadIdx.y,threadIdx.x)=c;

}

但是这种实现方式依旧存在访存的问题：由于每个线程还是计算block中的一个元素，还是存在对共享内存的反复读取，虽然共享内存的速度要比全局内存快很多，但是如果矩阵规模很大，那这一块的访存开销还是很大。

减少访存策略2：Thread分块

参考Block分块的思想，我们让每个线程也计算一个子块，Thread分块的基本思想与Block分块的思想相同：每个线程计算一个子块，计算的时候先将每个线程对应的A矩阵和B矩阵的共享内存数据读取到寄存器中，然后再计算，由于寄存器数量非常有限，所以这里也不能一次性加载到寄存器中，所以也需要分批加载，为了只读取一遍共享内存，我们需要采用矩阵外积的计算形式：也就是每个线程计算子块的时候，读取A矩阵的一列和B矩阵的一行到寄存器，然后计算外积，这样就只需要读取一遍共享内存就可以了。
在这里插入图片描述
代码实现如下：

// 分块参数
#define BM 128 // block子块大小
#define BN 128
#define BK 8
#define TM 8 // thread子块大小
#define TN 8

__global__ void BlockGEMM_V2(Matrix<float> A,Matrix<float> B,Matrix<float> C)
{
    // 每个线程的计算结果
    float c[TM][TN]={0.0};
    float a[TM]={0.0};
    float b[TN]={0.0};

    // 沿着K维度循环加载一个block中对应的A和B的数据到共享内存
    for(int i=0;i<A.cols/BK;++i)
    {
        // 每个block对应的全局内存中的A,B子块，即创建全局内存中A,B的view
        Matrix<float> ASub(A.data+blockIdx.y*BM*A.strideOfRow+i*BK,BM,BK,A.strideOfRow,A.strideOfCol);
        Matrix<float> BSub(B.data+i*BK*B.strideOfRow+blockIdx.x*BN,BK,BN,B.strideOfRow,B.strideOfCol);

        // 将Asub,BSub加载到共享内存
        // 以block为128，thread为8为例：由于一个block有16x16=256个线程，而ASub和BSub中一共有1024个元素，所以每个线程加载4个元素
        // 注意：这里需要将一维逻辑索引转换为多维逻辑索引：stardIndex->(stardIndex/cols,stardIndex%cols)
        __shared__ float A_Shared[BM][BK];
        __shared__ float B_Shared[BK][BN];
        int numberOfElementsPerThread=(BK*BM)/(blockDim.x*blockDim.y);// 每个线程需要读取多少数据
        int stardIndex=numberOfElementsPerThread*(threadIdx.y*blockDim.x+threadIdx.x);// stardIndex为每个线程读取的起始索引
        for(int threadIndex=0;threadIndex<numberOfElementsPerThread;++threadIndex)
        {
            int logicalIndex=stardIndex+threadIndex;
            A_Shared[logicalIndex/BK][logicalIndex%BK]=ASub(logicalIndex/BK,logicalIndex%BK);
            B_Shared[logicalIndex/BN][logicalIndex%BN]=BSub(logicalIndex/BN,logicalIndex%BN);
        }
        __syncthreads();

        // 每个thread对应的共享内存中的A_Shared,B_Shared的子块，即创建A_Shared,B_Shared的view
        Matrix<float> ASub_Shared((float *)A_Shared+threadIdx.y*TM*BK,TM,BK,BK,1);// 每个线程对应的共享内存中A和B的子块
        Matrix<float> BSub_Shared((float *)B_Shared+threadIdx.x*TN,BK,TN,BN,1);

        // 每个线程执行计算
        for(int k=0;k<BK;++k)
        {
            // 先将A的一列和B的一行加载到寄存器
            for(int m=0;m<TM;++m)
            {
                a[m]=ASub_Shared(m,k);
            }
            for(int n=0;n<TN;++n)
            {
                b[n]=BSub_Shared(k,n);
            }

            // 使用寄存器计算
            for(int m=0;m<TM;++m)
            {
                for(int n=0;n<TN;++n)
                {
                    c[m][n]+=a[m]*b[n];
                }
            }
        }
        __syncthreads();

    }

    // 将每个线程计算好的结果写回到C矩阵
    // CSub为每个线程对应的全局内存的C矩阵子块，创建C矩阵的view
    Matrix<float> CSub(C.data+((blockIdx.y*BM+threadIdx.y*TM)*C.strideOfRow+blockIdx.x*BN+threadIdx.x*TN),TM,TN,C.strideOfRow,C.strideOfCol);
    for(int m=0;m<TM;++m)
    {
        for(int n=0;n<TN;++n)
        {
            CSub(m,n)=c[m][n];
        }
    }

}

数据预取

通过前面的优化相对于NaiveGEMM已经有了很大提高了。下面我们分析一下BlockGEMM_V2的实现，在CUTLASS官方博客中的Software Pipelining一节中提到了BlockGEMM_V2这种实现方式存在这样的问题：每个线程按照“访存1—计算1—访存2—计算2—…—访存n—计算n”的顺序执行，这种执行方式每次计算单元都需要等待访存，而访存的延迟通常都是比较大的，所以这种实现会存在较大的访存延迟，为了减少访存延迟，Software Pipelining一节中提到了将下一次访存和上一次计算并行，这样可以掩盖访存的延迟。
在这里插入图片描述
这种优化方式也叫数据预取。代码实现如下：

// 分块参数
#define BM 128 // block子块大小
#define BN 128
#define BK 8
#define TM 8 // thread子块大小
#define TN 8

__global__ void BlockGEMM_V3(Matrix<float> A,Matrix<float> B,Matrix<float> C)
{
    // 每个线程的计算结果
    float c[TM][TN]={0.0};
    float a[TM]={0.0};
    float b[TN]={0.0};

    // 此时需要的共享内存是原来的2倍
    // 注意：读取和写入的时候第一个维度的索引是交错进行的
    __shared__ float A_Shared[2][BM][BK];
    __shared__ float B_Shared[2][BK][BN];

    // 预取(先读取第一个BK)
    Matrix<float> ASub(A.data+blockIdx.y*BM*A.strideOfRow+0*BK,BM,BK,A.strideOfRow,A.strideOfCol);
    Matrix<float> BSub(B.data+0*BK*B.strideOfRow+blockIdx.x*BN,BK,BN,B.strideOfRow,B.strideOfCol);
    int numberOfElementsPerThread=(BK*BM)/(blockDim.x*blockDim.y);
    int stardIndex=numberOfElementsPerThread*(threadIdx.y*blockDim.x+threadIdx.x);// stardIndex为每个线程读取的起始索引
    for(int threadIndex=0;threadIndex<numberOfElementsPerThread;++threadIndex)
    {
        int logicalIndex=stardIndex+threadIndex;
        A_Shared[0][logicalIndex/BK][logicalIndex%BK]=ASub(logicalIndex/BK,logicalIndex%BK);
        B_Shared[0][logicalIndex/BN][logicalIndex%BN]=BSub(logicalIndex/BN,logicalIndex%BN);
    }
    __syncthreads();

    // 沿着K维度循环加载剩下的数据
    int indexOfRead,indexOfWrite;
    bool indexFlag=false;// 辅助变量，用来计算索引
    for(int i=1;i<A.cols/BK;++i)
    {
        // 计算索引，indexOfRead和indexOfWrite每次循环会交替变换，i=1时为indexOfRead=0,indexOfWrite=1，i=2时为indexOfRead=1,indexOfWrite=0
        indexOfRead = (int)indexFlag; // 读索引，即本次循环读取A_Shared[indexOfRead,:,:]和B_Shared[indexOfRead,:,:]中的数据执行计算
        indexOfWrite = 1-indexOfRead; // 写索引，即预取下一次计算需要的数据到A_Shared[indexOfWrite,:,:]和B_Shared[indexOfWrite,:,:]中

        // 每个线程执行计算
        Matrix<float> ASub_Shared(((float *)A_Shared+indexOfRead*BM*BK)+threadIdx.y*TM*BK,TM,BK,BK,1);// 每个线程对应的共享内存中A和B的子块
        Matrix<float> BSub_Shared(((float *)B_Shared+indexOfRead*BK*BN)+threadIdx.x*TN,BK,TN,BN,1);
        for(int k=0;k<BK;++k)
        {
            // 先将A的一列和B的一行加载到寄存器
            for(int m=0;m<TM;++m)
            {
                a[m]=ASub_Shared(m,k);
            }
            for(int n=0;n<TN;++n)
            {
                b[n]=BSub_Shared(k,n);
            }

            // 使用寄存器计算
            for(int m=0;m<TM;++m)
            {
                for(int n=0;n<TN;++n)
                {
                    c[m][n]+=a[m]*b[n];
                }
            }
        }

        // 预取下个循环的数据
        Matrix<float> ASub(A.data+blockIdx.y*BM*A.strideOfRow+i*BK,BM,BK,A.strideOfRow,A.strideOfCol);
        Matrix<float> BSub(B.data+i*BK*B.strideOfRow+blockIdx.x*BN,BK,BN,B.strideOfRow,B.strideOfCol);
        int numberOfElementsPerThread=(BK*BM)/(blockDim.x*blockDim.y);
        int stardIndex=numberOfElementsPerThread*(threadIdx.y*blockDim.x+threadIdx.x);// stardIndex为每个线程读取的起始索引
        for(int threadIndex=0;threadIndex<numberOfElementsPerThread;++threadIndex)
        {
            int logicalIndex=stardIndex+threadIndex;
            A_Shared[indexOfWrite][logicalIndex/BK][logicalIndex%BK]=ASub(logicalIndex/BK,logicalIndex%BK);
            B_Shared[indexOfWrite][logicalIndex/BN][logicalIndex%BN]=BSub(logicalIndex/BN,logicalIndex%BN);
        }
        __syncthreads();

        // 设置flag
        indexFlag=!indexFlag;
    }

    // 计算最后一个BK
    {
        Matrix<float> ASub_Shared(((float *)A_Shared+indexOfWrite*BM*BK)+threadIdx.y*TM*BK,TM,BK,BK,1);// 每个线程对应的共享内存中A和B的子块
        Matrix<float> BSub_Shared(((float *)B_Shared+indexOfWrite*BK*BN)+threadIdx.x*TN,BK,TN,BN,1);
        for(int k=0;k<BK;++k)
        {
            // 先将A的一列和B的一行加载到寄存器
            for(int m=0;m<TM;++m)
            {
                a[m]=ASub_Shared(m,k);
            }
            for(int n=0;n<TN;++n)
            {
                b[n]=BSub_Shared(k,n);
            }

            // 使用寄存器计算
            for(int m=0;m<TM;++m)
            {
                for(int n=0;n<TN;++n)
                {
                    c[m][n]+=a[m]*b[n];
                }
            }
        }
    }

    // 将每个线程计算好的结果写回到C矩阵
    // CSub为每个线程对应的全局内存的C矩阵子块，创建C矩阵的view
    Matrix<float> CSub(C.data+((blockIdx.y*BM+threadIdx.y*TM)*C.strideOfRow+blockIdx.x*BN+threadIdx.x*TN),TM,TN,C.strideOfRow,C.strideOfCol);
    for(int m=0;m<TM;++m)
    {
        for(int n=0;n<TN;++n)
        {
            CSub(m,n)=c[m][n];
        }
    }

}

GEMM优化总结

我们用CUTLASS博客中的一张图总结一下GEMM优化的思路（即GEMM层次结构）：block分块->thread分块
在这里插入图片描述

性能测试

下面我们对比一下NaiveGEMM，BlockGEMM_V1，BlockGEMM_V2，BlockGEMM_V3和cublas的性能。

  1. 测试环境：P100，CUDA10.0
  2. 分块大小：BM=BM=128，BK=TM=TN=8
  3. 单位:ms
  4. 测试的时候，不考虑边界的处理，即M,N,K都是分块大小的整数倍，实验中A和B都设置为方阵，也就是M=N=K。
  5. 测试的时候BlockGemm_V1在128~1024之间选取的block大小为16x16，2048和4096选取的block大小为32x32
128256512102420484096
NaiveGemm0.030.080.504.2545.36402.47
BlockGemm_V10.020.060.251.729.4374.20
BlockGemm_V20.090.150.270.813.7726.51
BlockGemm_V30.090.150.270.773.6525.59
cublas0.020.020.050.352.0815.52

测试结果分析：

  1. 当矩阵规模较小的时候（<512）NaiveGemm的性能还不错，但是当矩阵规模逐渐变大之后，NaiveGemm性能就会显著下降
  2. 当矩阵规模较大的时候（>512），BlockGemm_V1相比于NaiveGemm性能有显著提高，说明优化全局内存的访存次数可以显著提高性能，BlockGemm_V2相比于BlockGemm_V1也有较大性能提升，说明优化共享内存的访存次数也可以显著提高性能，但是相比而言没有V1的提升幅度大
  3. 当矩阵规模较小的时候，比如128或者256，V2的性能是要比V1差的，由于V2每个block计算128x128的分块，所以在矩阵规模为128或者256的时候，参与计算的block数量很少，影响了并行性，虽然V2做了访存的优化，但是由于矩阵规模很小，这一块的优化提升并不明显，而V1由于每个block值计算16x16的分块，所以相比于V2，参与计算的block数量较多，所以性能要比V2好，对于V2如果修改分块参数为BM=BM,BK=16，TM=TN=1，则性能就与V1一致了
  4. 在P100，CUDA10.0环境中V3和V2性能差距很小，说明数据预取对性能提升很有限
  5. 从上面的性能数据来看，性能优化的主要收益还是来自于访存的优化。
  6. 本文中的GEMM优化与cublas的差距依旧很大！

最后总结一下，其实GEMM优化的核心还是在于访存的优化，即减少对延迟高的内存的重复访问，将需要重复访问的数据从延迟高的读取到延迟低的内存中执行计算，本文中GEMM的优化采用的是分块的策略实现了访存的优化。

完整测试代码见github项目: GEMM，欢迎star

结束语

由于自己也是刚接触GEMM的优化，所以对GEMM的优化研究并不是很深入，本文只是实现了一些最基本的GEMM优化方法，其实还可以做很多其他优化，比如bank冲突优化，而且本文并没有考虑到边界的处理。文中有什么不对的地方，欢迎批评指正。

2022-10-17 16:16:38
Last Updated: 2022-11-4 20:21:42

qianqing13579
关注
  • 3
    点赞
  • 踩
  • 26
    收藏
    觉得还不错? 一键收藏
  • 3
    评论
专栏目录
MegEngine| CUDA 矩阵乘法终极优化
旷视开发者技术社区
09-24 1166
前言单精度矩阵乘法（SGEMM）几乎是每一位学习 CUDA 的同学绕不开的案例，这个经典的计算密集型案例可以很好地展示 GPU 编程中常用的优化技巧，而能否写出高效率的 SGEMM Ker...
dev c++无法识别汇编代码_新手福利丨超详细的Tengine GEMM矩阵乘法汇编教程
weixin_39540744的博客
12-15 413
很多刚入门Tengine的开发者想研读Tengine汇编代码，却苦于没有好的汇编入门教程，没有大神带入门，自己看又看不懂，怎么办？福利来了，Tengine带来了一份超详细的gemm汇编教程。GEMM简介什么是GEMM? 它的英文全称是 GEneral Matrix to Matrix Multiplication (通用矩阵的矩阵乘法)，Gemm在神经网络的计算中占据很重要的位置。Why...
CUDA算子优化：矩阵乘GEMM优化（一）
qq_56047026的博客
06-12 836
采用最navie的GEMM算法，在GPU中，一共开启m*n个线程，每个线程需要读取矩阵A的一行和矩阵B的一列，而后将计算结果写回矩阵C中。其中，M=m/bm，N=n/bn，K=k/bk。当进行K轮迭代中某一轮迭代时，GPU将维度为bm*bk，bk*bn的小矩阵块存储到shared memory中，而后各个线程将shared memory中的数据存入register中进行计算。在完成以上的过程后，对于一个线程而言，它现在拥有：rm个A矩阵的寄存器值，rn个B矩阵的寄存器值，一级rm*rn个C矩阵的寄存器值。
CUDA算子优化：矩阵乘GEMM优化(二)
最新发布
qq_56047026的博客
06-13 961
A_TILE_ROW_START代表在这个维度为bm*bk的数据块中，当前线程需要搬运的数据的竖向坐标，而A_TILE_COL代表需要搬运的数据的横向坐标。那么我们需要开启（2048/128）*（2048/128）=256个block，每个block里面有（128/8）*（128/8）=256个线程，每个线程需要负责计算C矩阵中8*8=64个元素的结果，每个block负责256*64=16384个元素的结果。在一轮迭代中需要使用bm*bk的数据，为了加快后续的访存，所以需要进行一次转置。
GEMM算法及优化流程详解
热门推荐
独行侠
02-16 2万+
前言 神经网络前向耗时主要由卷积的耗时决定，参考賈杨青毕业论文，那么如何对卷积加速便成了重要的一个点，主流的加速方法有 以下几种： im2col+GEMM：目前几乎所有的主流计算框架包括 Caffe, MXNet 等都实现了该方法. 该方法把整个卷积过程转化成了GEMM过程，而GEMM在各种 BLAS 库中都是被极致优化的，一般来说，速度较快。 Winograd： Winograd 是存在...
[CUDA 学习笔记] GEMM 优化: 双缓冲 (Prefetch) 和 Bank Conflict 解决
LostUnravel的博客
04-29 1760
本文主要是对以及三篇文章相关内容的学习整理.本文可以作为这篇文章的补充. 可以看作文中的后续优化. 在先前的这篇文章中, 作者原文中并没有针对 GEMM 的"双缓冲"和"Bank Conflict"进行具体讲解, 虽然译者根据原文作者提供的代码进行了这两部分的补充, 但译文中也提到了代码中存在的一些问题. 上述提到的三篇文章恰好提到了这部分内容, 通过比较, 笔者认为其提供了更好的解决方案, 值得学习.
令人头秃的cudaTensorCoreGemm详解
jianjianheng的博客
08-27 5226
我曾被cudaTensorCoreGemm的代码苦苦折磨了一个月，官方文档也不够详细，代码虽然很多注释也不太好理解，于是决定自己写一个cudaTensorCore的详解。本文主要对cudaTensorCore中的compute_gemm函数进行分析。 Compute_gemm函数流程图如下图所示： Cuda wmma API主要用于使用SM中WARP和TensorCore加速矩阵乘加运算，标准公式如下（大写字母均为矩阵）： D=α×A×B+β×C 其中A矩阵为M×K，B...
CUDA矩阵乘法GEMM优化：全局内存-共享内存-寄存器优化，以及数据预存取优化
weixin_55035144的博客
06-21 1004
可以将 A 和 B 矩阵先搬运到 Shared Memory（SM 中低延迟的 on-chip memory，block 内线程共享，附 NVIDIA GPU 内存结构图）中降低访存的开销，这的确是一个很好的思路，但是这只能将访存代价从几百 cycle 降低到几十 cycle，并不改变问题的本质。相对于之前的情况，计算指令的占比大大提高了。问题的关键在于主体循环由两条 Load 指令与一条 FMA 指令构成，计算指令只占总体的 1/3，计算访存比过低，最终导致了访存延迟不能被隐藏，从而性能不理想。
cuda-使用cuda并行加速实现之gemv.zip
04-24
2. CUDA编程模型：理解流式多处理器（SM）、线程块和线程的概念。 3. 内存管理：CUDA内存层次（全局、共享、常量、纹理、对齐），以及cudaMalloc、cudaFree、cudaMemcpy等函数的使用。 4. 内核函数：编写和执行CUDA...
用CUDA（显卡）解线性方程组
04-05
首先，我们需要理解CUDA编程的基本概念，如设备内存模型、线程块和网格、同步机制以及数据传输等。例如，线程在GPU上是按照线程块和网格组织的，每个线程块内的线程可以高效地共享内存，而不同线程块间的通信则通过...
cuda_9.0.176.2补丁就是他.zip
10-20
CUDA是NVIDIA公司推出的一种并行计算平台和编程模型，主要应用于高性能计算和图形处理。...对于开发者来说，熟悉CUDA编程模型和错误调试方法是至关重要的，这样才能有效利用GPU的计算能力，避免或快速解决类似的问题。
CUDA矩阵操作
12-04
CUDA矩阵操作是基于...通过理解和应用CUDA编程模型，开发者可以构建出比传统CPU实现更快、更节能的解决方案，尤其在处理大规模矩阵问题时。对于学习和使用CUDA的人来说，理解这些基本概念和优化策略是至关重要的。
cublas_benchmarks：一些测试gemm和gemv性能的基准
02-17
在CUDA编程环境中，CuBLAS（CUDA BLAS）是一个由NVIDIA提供的库，它利用GPU的强大计算能力来加速线性代数操作，如通用矩阵乘法（GEMM）和通用向量乘法（GEMV）。这些操作是许多科学计算、机器学习和深度学习算法的...
通用矩阵乘（GEMM）优化与卷积计算
u013250861的博客
11-26 1433
前两节探讨的是传统优化矩阵乘的两者思想，近年来兴起的神经网络中广泛应用的矩阵乘也可用类似方法优化（实际上深度学习软件可以直接使用已有的数学加速库，如 BLAS）。随着技术的演进，神经网络技术出现了一个重要的方向——神经网络量化。在常规的数据规模中，由于左侧会发生太多的高速缓存缺失，又由于矩阵乘这样的计算对数据的访问具有很高的重复性，将它重排成右侧的内存布局减少高速缓存缺失，可显著地改进性能。到目前为止（2019 年中），QNNPACK 是已公开的，用于移动端（手机）的，性能最优的量化神经网络加速库。
BLAS之GEMM和GEMV
永远飞翔的鸟
10-17 7844
BLAS简介 BLAS全称是Basic Linear Algebra Subprograms是规定了一套低级的执行常见线性代数操作的规范。其实现经常针对特殊的机器进行优化，比较著名的·BLAS库有ACML, ATLAS, MKL, OpenBLAS。许多常见的数值软件均采用兼容BLAS规范的实现库来进行线性代数计算，比如Matlab, Numpy, Mathematica`。 其中，Level 1 BLAS主要提供向量操作 Level 2 BLAS提供矩阵向量操作（gemv) gemv
CUDA基础（四）杂。GEMM优化思想、内存事务。
IsayIwant的博客
07-27 346
为操作码（Opcode，当然它本质上还是指的FFMA对应的编码）, 后面的R2, R4, R2, R5表示参与操作的通用寄存器（General Purpose Register，简称GPR，有时也直接叫Register），称为操作数（Operand）。每一个内存事务，简单说就是一个邮递员，他每次可以携带一封信从A到B，信最大可以携带128字节的信息量，通勤一次时间200天，所以你是愿意写一封128字节的信，还是写32封4字节的信，方案不言而喻。所以说，GEMM优化的重要性，怎么突出都不过分。
CUDA 并行计算优化策略总结
qq_29462849的博客
12-29 877
作者 | LustofLife@知乎来源 | https://zhuanlan.zhihu.com/p/297201517编辑 | 极市平台导读并行计算为了提高算法运行效率，本文通过以矩...
GEMM矩阵相乘与深度学习
weixin_43728590的博客
07-08 4248
GEMM矩阵相乘与深度学习1. GEMM矩阵相乘1.1 GEMM算法基础1.2 GEMM算法优化1.2.1 循环重排充分利用缓存1.2.2 平铺（Tiling）充分利用缓存1.2.3 展开（Unrolling）1.2.4 内存对齐1.2.5 向量化1.2.6 矩阵分块1.2.7 双缓冲2. GEMM与卷积计算2.1 计算卷积的方法2.2 Img2col2.2.1 CNN中张量的存储2.2.2 卷积运算转化为GEMM 1. GEMM矩阵相乘 1.1 GEMM算法基础   GEMM（General Matrix
Blas GEMM launch failed
04-03
As an AI language model, I do not have access to your system or the context of your question, so I can only provide a general answer. Blas GEMM stands for Basic Linear Algebra Subprograms General Matrix-Matrix multiplication, which is a common operation in linear algebra and machine learning. If you are encountering a "Blas GEMM launch failed" error, it could mean that there is an issue with your system's BLAS library or that the input matrices are not compatible for the operation. To troubleshoot the issue, you can try the following: 1. Check that your system has a BLAS library installed and that it is up to date. 2. Verify that the input matrices have the correct dimensions for matrix multiplication, i.e., the number of columns in the first matrix must match the number of rows in the second matrix. 3. Check for any memory issues or conflicts that may be preventing the operation from launching. 4. Consult the documentation or support forums for the specific software or library you are using for additional guidance on troubleshooting the error. In general, the "Blas GEMM launch failed" error can be caused by a variety of factors, so it may require some experimentation and investigation to pinpoint the exact issue.
评论 3
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值