CUDA编程之GEMM优化

置顶 已于 2022-11-04 20:22:04 修改
阅读量2.4k 收藏 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 1145
前言单精度矩阵乘法（SGEMM）几乎是每一位学习 CUDA 的同学绕不开的案例，这个经典的计算密集型案例可以很好地展示 GPU 编程中常用的优化技巧，而能否写出高效率的 SGEMM Ker...
如何在 CPU 上优化 GEMM
吴建明wujianming_110117
10-05 362
如何在 CPU 上优化 GEMM (TL;DR) TVM 提供抽象接口，允许用户分别描述算法和算法的实施组织（所谓的调度）。通常，在高性能调度中编写算法，会破坏算法的可读性和模块化。尝试各种看似有前途的调度也很耗时。在 TVM 的帮助下，可以有效地尝试这些调度，提高性能。 将演示如何使用 TVM 优化矩阵乘法，通过简单地添加 18 行额外代码，实现比基线快 200 倍。 在 CPU 上执行的密集计算应用程序，有两个重要的优化： 提高内存访问的缓存命中率。复杂的数值计算和热点内存访问，都可以通过高缓存命中率
CUDA算子优化：矩阵乘GEMM优化（一）
qq_56047026的博客
06-12 794
采用最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 优化: 双缓冲 (Prefetch) 和 Bank Conflict 解决
LostUnravel的博客
04-29 1464
本文主要是对以及三篇文章相关内容的学习整理.本文可以作为这篇文章的补充. 可以看作文中的后续优化. 在先前的这篇文章中, 作者原文中并没有针对 GEMM 的"双缓冲"和"Bank Conflict"进行具体讲解, 虽然译者根据原文作者提供的代码进行了这两部分的补充, 但译文中也提到了代码中存在的一些问题. 上述提到的三篇文章恰好提到了这部分内容, 通过比较, 笔者认为其提供了更好的解决方案, 值得学习.
CUDA算子优化：矩阵乘GEMM优化(二)
最新发布
qq_56047026的博客
06-13 892
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的数据，为了加快后续的访存，所以需要进行一次转置。
CUDA 并行计算优化策略总结
qq_29462849的博客
12-29 870
作者 | LustofLife@知乎来源 | https://zhuanlan.zhihu.com/p/297201517编辑 | 极市平台导读并行计算为了提高算法运行效率，本文通过以矩...
GEMM实现（NEON+CUDA）
qq_45642593的博客
12-13 1564
在Fortran编程当中数组是默认按列存储的，而在C/C++当中通常是按行存储的。A的第一行先和B的第一行按照A的通道为1的元素来进行广播相乘，因为B0就是矩阵B的第一行，按照矩阵的乘法，只能和A的第一行的第一个元素相乘，然后累加到寄存器c0中，这样将B的四个行遍历结束后，就能计算出C矩阵中的一行了。然后将上面的思想用到一张640*640的数组中计算，其实就是对矩阵A的行进行4个为一组的分行块，对矩阵B的列进行分列块，每一次分好的块进行4*4的矩阵计算，就能够迭代出整个矩阵C的结果了。
产线电池电压ADC校准方案
爱吃肉的猫
12-23 1617
校准方法如下： 外部给固定电压V11，设备ADC采到的值为V1；外部再给出固定电压V22，设备内部ADC采到的值为V2；则： V11 = Gain * V1 + Offset； V22 = Gain * V2 + Offset； 利用以上两个计算式，计算出Gain和 Offset； V矫正后 = Gain * V采集值 + Offset； ...
用CUDA（显卡）解线性方程组
04-05
首先，我们需要理解CUDA编程的基本概念，如设备内存模型、线程块和网格、同步机制以及数据传输等。例如，线程在GPU上是按照线程块和网格组织的，每个线程块内的线程可以高效地共享内存，而不同线程块间的通信则通过...
cuda_9.0.176.2补丁就是他.zip
10-20
CUDA是NVIDIA公司推出的一种并行计算平台和编程模型，主要应用于高性能计算和图形处理。...对于开发者来说，熟悉CUDA编程模型和错误调试方法是至关重要的，这样才能有效利用GPU的计算能力，避免或快速解决类似的问题。
CUDA矩阵操作
12-04
CUDA矩阵操作是基于...通过理解和应用CUDA编程模型，开发者可以构建出比传统CPU实现更快、更节能的解决方案，尤其在处理大规模矩阵问题时。对于学习和使用CUDA的人来说，理解这些基本概念和优化策略是至关重要的。
拉伸：为GEMM和张量收缩扩展GPU性能
02-23
GEMM，全称Generalized Matrix Multiply，是线性代数中最基础的操作之一，广泛应用于机器学习、深度神经网络(DNN)的卷积操作以及其他高性能计算领域。 GEMM不仅限于基本的二维矩阵乘法，还涵盖了批量矩阵乘法...
cublas_benchmarks：一些测试gemm和gemv性能的基准
02-17
在CUDA编程环境中，CuBLAS（CUDA BLAS）是一个由NVIDIA提供的库，它利用GPU的强大计算能力来加速线性代数操作，如通用矩阵乘法（GEMM）和通用向量乘法（GEMV）。这些操作是许多科学计算、机器学习和深度学习算法的...
【架构分析】GPU执行GEMM矩阵运算实例演示
HaoBBNuanMM的博客
06-09 444
Cutlass是 NVIDIA 提供的一套用于高效实现矩阵乘法和卷积操作的 C++ 库。它以 CUDA 为基础，提供了高度优化的数学运算，尤其适用于GPU上的高性能并行计算。本文以GEMM矩阵运算作为实例，展示Cutlass在GPU上执行GEMM运算的过程。
YOLOv3目标检测：原理与源码解析
06-07
Linux创始人LinusTorvalds有一句名言：Talk is cheap, Show me the code.（冗谈不够，放码过来！）。 代码阅读是从入门到提高的必由之路。尤其对深度学习，许多框架隐藏了神经网络底层的实现，只能在上层调包使用，对其内部原理很难认识清晰，不利于进一步优化和创新。  YOLOv3是一种基于深度学习的端到端实时目标检测方法，以速度快见长。YOLOv3的实现Darknet是使用C语言开发的轻型开源深度学习框架，依赖少，可移植性好，可以作为很好的代码阅读案例，让我们深入探究其实现原理。  本课程将解析YOLOv3的实现原理和源码，具体内容包括： YOLO目标检测原理  神经网络及Darknet的C语言实现，尤其是反向传播的梯度求解和误差计算 代码阅读工具及方法 深度学习计算的利器：BLAS和GEMM GPU的CUDA编程方法及在Darknet的应用 YOLOv3的程序流程及各层的源码解析本课程将提供注释后的Darknet的源码程序文件。  除本课程《YOLOv3目标检测：原理与源码解析》外，本人推出了有关YOLOv3目标检测的系列课程，包括：   《YOLOv3目标检测实战：训练自己的数据集》  《YOLOv3目标检测实战：交通标志识别》  《YOLOv3目标检测：原理与源码解析》  《YOLOv3目标检测：网络模型改进方法》 建议先学习课程《YOLOv3目标检测实战：训练自己的数据集》或课程《YOLOv3目标检测实战：交通标志识别》，对YOLOv3的使用方法了解以后再学习本课程。
cuda三维索引与一维索引（一）
qq_43594926的博客
03-22 1028
三维网格中每个网格的坐标为 (i,j,k)，其中 i 表示沿 x 轴的索引，j 表示沿 y 轴的索引，k 表示沿 z 轴的索引。nx、ny、nz 分别为数组在 x、y、z 轴上的长度其中，blockIdx.x、blockDim.x、threadIdx.x 表示当前线程所在的块在 x 轴上的索引、每个块在 x 轴上的线程数、当前线程在所在块中在 x 轴上的索引，其他变量类似。这些变量通常在 CUDA 编程中使用。则 (i, j, k) 在一维数组中的索引为。
【BBuf的CUDA笔记】七，总结 FasterTransformer Decoder(GPT) 的cuda相关优化技巧
I good vegetable a!
02-05 1808
解读：https://github.com/NVIDIA/FasterTransformer/blob/main/docs/gpt_guide.mdGPT 是 Decooding 模型的一种变体，没有 Encoder 模块，没有交叉多头注意力模块，使用 GeLU 作为激活函数。2020 年，OpenAI 在他们的论文中表明，使用非常庞大的模型和大量的训练数据可以显著提高 GPT 模型的容量。但是，不可能将这样的模型放入单个 GPU 中。
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、付费专栏及课程。

余额充值