CUDA on Platform 学习笔记4--存储单元及矩阵乘法

CV温故知新

已于 2022-07-08 16:00:27 修改

阅读量322

点赞数

分类专栏： CUDA教程

于 2022-07-08 11:29:54 首次发布

本文链接：https://blog.csdn.net/lansebingxuan/article/details/125662153

版权

深度学习

CUDA教程专栏收录该内容

6 篇文章 2 订阅

订阅专栏

本文是参加2022 CUDA on Platform 线上训练营学习笔记，感谢NVIDIA各位老师的精彩讲解!

1. GPU的存储单元

GPU的存储单元分为两大类：

板子上芯片周围的显存颗粒(on board)，读取速度相对慢，如下图中的local memory，global memory，constant memory，texture memory。
在GPU芯片内部(on chip)，读取速度相对快，如下图中的registers，shared memory。

下图中的箭头双向表示可以读写，单向表示只能读。这些memory可以进一步细分：
1）R/W可读可写memory：

registers、local memory：线程私有memory，每个线程私有访问。
shared memory：一个block内的线程都可以访问，可以数据共享通信。
global memory：每个线程都可以读写。

2）R 只读memory：constant memory，texture memory是每个线程都可以读。
global memory，constant memory，texture memory和主机之间都可以通信读写，通常显卡说明书写的显存大小，就是global memory
在这里插入图片描述

2. GPU存储单元的分配与释放

CPU和GPU之间的memory copy函数
cudaMemcpy(void *dst，const void *src，size_t count，cudaMemcpyKind kind），参数含义为：
1）dst: destination memory address
2）src: source memory address
3） count: size in bytes to copy
4）kind: direction of the copy
拷贝方向参数cudaMemcpyKind的类型有下面4种：
1）cudaMemcpyHostToDevice
2）cudaMemcpyDeviceToHost
3）cudaMemcpyDeviceToDevice
4）cudaMemcpyHostToHost
申请GPU存储单元
当我们要为一个方阵M(m * m)申请GPU的存储单元时，使用下面的函数：
cudaMalloc((void**) &d_m，sizeof(int) * m * m)，参数含义如下：
1）d_m：指向存储在Device端数据的地址的指针
2）sizeof(int) * m * m：存储在Device端空间的大小
这里需要注意(void)双重指针**，由于cudaMalloc函数返回类型是错误类型cudaError_t，因此分配的内存地址只能由参数传入修改并返回。当在函数里修改这个值的时候需要传入这个值的地址，因此需要一个指针指向这个地址，而分配内存时候返回一个指针，指向分配地址，因此该参数为指向指针的指针。
释放GPU申请的存储单元函数
cudaFree(d_m);
d_m：指向存储在Device端数据的地址的指针
从CPU内存传输到GPU存储单元
cudaMemcpy(d_m，h_m，sizeof(int）* m * m，cudaMemcpyHostToDevice），各参数的设置为：
1）d_m：传输的目的地，GPU存储单元
2）h_m：数据的源地址，CPU存储单元
3）sizeof(int）* m * m：数据传输的大小
4）cudaMemcpyHostToDevice：数据传输的方向，CPU到GPU
从GPU存储单元传输到CPU内存
当我们完成计算后将数据从GPU存储单元传输到CPU的内存时函数如下：
cudaMemcpy(h_c, d_c, sizeof(int) * m * k，cudaMemcpyDeviceToHost)，注意这里的传输方向。

矩阵相乘样例

1. CPU计算方式
假如要计算矩阵相乘P = M* N ，其中M矩阵维度为m* n, N矩阵维度为 n* k，则P矩阵维度为m * k。CPU执行该程序代码如下：

void cpu_matrix_mult(int * h_m, int * h_n, int * h_result,int m, int n, int k)
{
    for (int i=0; i<m; ++i) 
    {
        for (int j=0; j<k; ++j)
        {
            int tmp = 0;
            for (int h=0; h<n; ++h)
            {
                tmp+=h_m[i * n + h] * h_n[h * k + j];
            }
            h_result[i * k + j]=tmp;
        }
    }
}

其中，i是M矩阵的行索引，j是N矩阵中的列索引，h是N矩阵的行索引， i和j同时也是P矩阵的元素索引，具体计算过程为：
1）M矩阵的第i行乘以N矩阵中的第j列，得到P矩阵的第i行j列个数据（n次乘加）
2）步骤1中N矩阵j从0到k-1重复k次，得到P矩阵的第i行k个列的所有数据（k次乘加）
3) M矩阵的i从0到m-1重复m次，得到P矩阵的m个行，k个列的所有数据（m次乘加）
因此，执行一次矩阵相乘需要计算m* n * k次乘加。
2. GPU优化方式
CUDA优化是用空间换时间，去掉步骤2）步骤3）两维循环，将矩阵每个值由二维表示转化为一维表示，每个位置在全局有唯一索引。
在这里插入图片描述
左图是一个grid，每个位置索引是(ix,iy)，其中ix=blockIdx.x* blockDim.x + threadIdx.x， iy=blockIdx.y* blockDim.y + threadIdx.y。右图中，一个grid中，x方向2个block，y方向3个block，共6个block，x方向8个线程，nx=8。每个block中，x方向4个线程，y方向2个线程，共8个线程。
对于红框中的元素，y方向block在grid中的索引是2，blockIdx.y=2；对于每个block，y方向2个线程，因此blockDim.y是2，1个block中，该位置的线程索引是1，threadIdx.y=1，x方向计算同理，
因此 iy=blockIdx.y* blockDim.y + threadIdx.y=(2* 2) +1=5; ix=blockIdx.x* blockDim.x + threadIdx.x = (14 + 0)=4,
grid中第iy行第ix列转换到一维索引=iy nx + ix=5 * 8+4=44。

下面是将一个grid用二维索引详细表示，第一张图坐标是相对于每个block的索引值
在这里插入图片描述
当把整个grid对应到一个矩阵的时候，它就像下面的样子，每个坐标是相对于整个grid的索引值。

上述转换和具体线程对应图如下：

一个kernel grid计算Pd
每个线程计算Pd的一个元素
每个线程
1）读入矩阵Md的一行
2）读入矩阵Nd的一列
3）为每对Md和Nd元素执行一次乘法和加法

下面是具体代码实现：

__global__void gpu_matrix_mult(int*M, int*N,int*P,int m_size,int n_size,int k_size)
{
    int row = threadIdx.y + blockIdx.y + blockDim.y;
    int col = threadIdx.x + blockIdx.x + blockDim.x;
    int sum=0;
    if(col < k_size && row < m_size)
    {
        for(int i = 0; i < n_size; i++)
        {
            sum += M[row * n_size + i] * N[col + i* k_size];
        }
        P[row * k_size + col] = sum;
    }
}