【并行计算】GPU,CUDA

发布时间:2023年12月30日

一、CUDA层次结构

1.kernel核函数

一个CUDA程序是一个kernel核函数被GPU的多个计算单元并行执行的过程,CUDA给了如下抽象

dim3 threadsPerBlock(4, 3, 1);
dim3 numBlocks(3, 2, 1);
matrixAdd<<<numBlocks, threadsPerBlock>>>(A, B, C);

?

2.Grid,Block,Thread

这样启动核函数,根据CUDA的抽象,就会有下面这样的运行模式,<<<>>>中间的两个参数numblocks和threadsPerBlock都是三维的变量,给予程序员设计的便利。

每个thread就是一个实际的核函数在运行,核函数可以根据当前的blockIdx,threadIdx来获得当前核函数所在的三维坐标位置。

int index = blockIdx.x * blockDim.x + threadIdx.x;

?

3.Streaming Multiprocessor(SM),warp

每个Block会分给一个SM(Streaming Multiprocessor),一个SM可以理解成一个有很多核的处理单元,并且有一个共享内存,下面看看一个SM内部如何工作。

下面这个图是一个典型的SM内部,每个黄方框都是一个SIMD单元,他们共享一个内存,左边的warp是实际分配给这些SIMD单元的任务,一个warp是一些线程的集合,CUDA用行优先的逻辑将一个block里的thread分配给warp,注意CUDA这里dim这个东西横纵坐标跟别的不太一样,如下图,他是Y是行号,X是列号。

在CUDA文档中,有讲到是根据线程id来顺序连续分配的,线程id计算方式如下

对于1维的来说,1维的x就是线程id

对于2维的来说,id是x + y Dx,y是行号,x是列号,所以就是行号乘一行的数量再加上列号。

对于3维的来说,id是x + y Dx + z Dx Dy,那就是高(z)乘上一个面的线程数,再加上y乘上行长在加上x。

所以总结来说,就是先分配面,然后在面上行优先分配。

一个warp通常是32个thread来执行SIMD指令,因为每个线程都是同样的核函数。但这里其实会有一个问题,那就是条件分支可能会不一样,最大的效率在这32个线程都执行相同的条件分支时达到,因为不同的分支会导致simd单元先执行一部分,而另一部分会等这部分执行完在执行。

所以一个warp才类似于操作系统中的一个线程,GPU会将warp视为线程来做硬件多线程调度。

看左边这一堆warp,存的就是每个warp的运行时状态,这里面包含了每个warp独立的寄存器、PC等东西,所以这里GPU做的硬件多线程就类似于一种超线程技术,使用多套上下文,使上下文切换没有开销。

二、CUDA内存层次结构

从最快的每个thread私有的内存,然后是整个块共享的一片内存,然后到整个GPU共享的全局内存。

一个值得注意的点,当一个warp访问内存中连续的地址时,会做块读取/写入,一次性将一个块内容读取/写入,所以如果让一个warp内的线程具有连续的内存访问模式,是比较好的,结合刚才的,如果也有同样的条件分支,那更好了。

三、一个矩阵乘法的优化例子

1.最基本的

直接A的行乘B的列相加,这会导致B的内存访问模式是跳跃的,不缓存友好。

2.预转置

那么就把B提前转置了,这样A和B都可以一行一行的访问了。

可以看到有一定的优化了

?

3.变成CUDA代码

最基础的版本,我们让C结果矩阵的每一个元素都用一个核函数来算结果,i和j就是C矩阵的i和j,我们直接将整个grid,映射成一个二维矩阵,那么横坐标i就是先拿块id的y乘上块的长度再加上块里面线程的横坐标y。纵坐标也类似。

___global__ void CUDASimpleKernel(int N, float *dmatA, float *dmatB, float *dmatC)
{
    int i = blockIdx.y * blockDim.y + threadIdx.y;
    int j = blockIdx.x * blockDim.x + threadIdx.x;
    if (i >= N || j >= N)
        return;
    float sum = 0.0;
    for (int k = 0; k < N; k++)
    {
        sum += dmatA[RM(i, k, N)] * dmatB[RM(k, j, N)];
    }
    dmatC[RM(i, j, N)] = sum;
}

然后i,j确定下来后,就去用k遍历A矩阵的一行和B矩阵的一列来计算结果元素。

当然,要变成CUDA代码还需要一些初始化的host代码。

首先要在GPU上分配内存,然后Memcpy过去

然后初始化块的数量和块的大小,就可以启动核函数了

然后算完之后再Memcpy回CPU

最后别忘了free掉GPU上用的内存

void CUDAMultMatrixSimple(int N, float *dmatA, float *dmatB, float *dmatC)
{
    dim3 threadsPerBlock(LBLK, LBLK);
    dim3 blocks(updiv(N, LBLK), updiv(N, LBLK));
    CUDASimpleKernel<<<blocks, threadsPerBlock>>>(N, dmatA, dmatB, dmatC);
}

void CUDAMultiply(int N, float *aData, float *bData, float *cData)
{
    float *aDevData, *bDevData, *cDevData;
    CUDAMalloc((void **)&aDevData, N * N * sizeof(float));
    CUDAMalloc((void **)&bDevData, N * N * sizeof(float));
    CUDAMalloc((void **)&cDevData, N * N * sizeof(float));
    CUDAMemcpy(aDevData, aData, N * N * sizeof(float), CUDAMemcpyHostToDevice);
    CUDAMemcpy(bDevData, bData, N * N * sizeof(float), CUDAMemcpyHostToDevice);

    CUDAMultMatrixSimple(N, aDevData, bDevData, cDevData);

    CUDAMemcpy(cData, cDevData, N * N * sizeof(float), CUDAMemcpyDeviceToHost);

    CUDAFree(aDevData);
    CUDAFree(bDevData);
    CUDAFree(cDevData);
}

好的,这有一个巨额的提升。

4. 考虑一个情况

刚才的i和j计算的代码变成这样,效果会变差十多倍。为什么呢

int i = blockIdx.x * blockDim.x + threadIdx.x;
int j = blockIdx.y * blockDim.y + threadIdx.y

想想内存访问模式的变化。刚才的代码,一个block里的一个warp是横着连续的,

横着连续说明他们的i一样,j连续,这说明,在对A矩阵的访问上,一直用的都是同一行,是内存中同一个连续的位置,可以进行块读。对B矩阵的访问上,是一列一列访问的,但是整个warp所需要访问的内存是连续的,所以也可以进行块读。

然后,对于写,是写C的连续的位置,因为是横着的,所以可以进行块写。

而新的代码

i是列号乘块的纵长,再加上块里的线程纵位置,也就是i和j对比刚才互换了,这样会导致什么,同一个warp里计算的是C矩阵纵向的元素。C矩阵纵向的元素,对于A,是不同的行,这样warp内整体也是连续的,可以进行块读,对于B,是同一列,这里读是不能块读的,因为内存是不连续的。

再看写,是竖着写的,所以写的也是C的不连续的位置,这样写也不能进行块写。

综上,这两个就差在一个块写和块读上了。

文章来源:https://blog.csdn.net/TheSysy/article/details/135302645
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。