🎉感谢来自NVIDIA企业级开发者社区的何琨(Ken)老师提供的资料和细致耐心的讲解
本篇不涉及存储单元的加速优化,重点演示最基础的内存申请释放、数据拷贝传输
存储单元
(这块并没有完全学完,Constant memory和Texture memory存储介绍内容来源于搜索引擎)
Register
寄存器是GPU最快的memory,寄存器变量是每个线程私有的,一旦thread执行结束,寄存器变量就会失效。寄存器是稀有资源。在Fermi上,每个thread限制最多拥有63个register,Kepler则是255个。让自己的kernel使用较少的register就能够允许更多的block驻留在SM中,也就增加了Occupancy,提升了性能。
寄存器不像共享内存、设备内存、常量内存和纹理内存显式声明,它们的声明不加任何限定符,就如普通变量一样,如 int a
。
寄存器主要承担Scalar variables(标量变量)和编译时已知的静态索引的数组。
Local memory
有时候,如果Register不够用了,那么就会使用Local memory来代替这部分寄存器空间。
除此外,下面几种情况,编译器可能会把变量放置在Local memory:
- 编译期无法决定确切值的本地数组。
- 较大的结构体或者数组,也就是那些可能会消耗大量Register的变量。
- 任何超过寄存器限制的变量。
在Local memory中的变量本质上跟Global memory在同一块存储区。所以,Local memory有很高的atency和较低的bandwidth。在CC2.0以上,GPU针对Local memory会有L1(per-SM)和L2(per-device)两级cache。
Shared memory
用__shared__修饰符修饰的变量存放在Shared memory。因为Shared memory是on-chip的,他相比Local Memory和Global memory来说,拥有高的多bandwidth和低很多的latency。他的使用和CPU的L1 cache非常类似,但是他是programmable的。
按惯例,像这类性能这么好的memory都是有限制的,Shared memory是以block为单位分配的。我们必须非常小心的使用Shared memory,否则会无意识的限制了active warp的数目。
不同于Register,Shared memory尽管在kernel里声明的,但是他的生命周期是伴随整个block,而不是单个thread。当该block执行完毕,他所拥有的资源就会被释放,重新分配给别的block。
Shared memory是thread交流的基本方式。同一个block中的thread通过Shared memory中的数据来相互合作。获取Shared memory的数据前必须先用__syncthreads()同步。
Global memory
Global Memory是空间最大,latency最高,GPU最基础的memory。“global”指明了其生命周期。任意SM都可以在整个程序的生命期中获取其状态。global中的变量既可以是静态也可以是动态声明。可以使用__device__修饰符来限定其属性。Global memory的分配就是之前频繁使用的cudaMalloc
,释放使用cudaFree
。Global memory驻留在Device memory,可以通过32-byte、64-byte或者128-byte三种格式传输。这些memory transaction必须是对齐的,也就是说首地址必须是32、64或者128的倍数。优化memory transaction对于性能提升至关重要。当warp执行memory load/store时,需要的transaction数量依赖于下面两个因素:
- Distribution of memory address across the thread of that warp 就是前文的连续
- Alignment of memory address per transaction 对齐
一般来说,所需求的transaction越多,潜在的不必要数据传输就越多,从而导致throughput efficiency降低。
Constant memory
Constant Memory驻留在Device Memory,并且使用专用的Constant cache(per-SM)。该Memory的声明应该以__connstant__修饰。Constant的范围是全局的,针对所有kernel,对于所有CC其大小都是64KB。在同一个编译单元,Constant对所有kernel可见。
当一个warp中所有thread都从同一个Memory地址读取数据时,Constant Memory表现最好。例如,计算公式中的系数。如果所有的thread从不同的地址读取数据,并且只读一次,那么Constant memory就不是很好的选择,因为一次读Constant memory操作会广播给所有thread知道。
Texture memory
Texture Memory驻留在Device memory中,并且使用一个只读cache(per-SM)。Texture Memory实际上也是Global Memory在一块,但是他有自己专有的只读cache。这个cache在浮点运算很有用。Texture Memory是针对2D空间局部性的优化策略,所以thread要获取2D数据就可以使用Texture Memory来达到很高的性能。
内存的分配与回收
CPU内存
- malloc()
- memset()
- free()
GPU内存
- cudaMalloc()
- cudaMemset()
- cudaFree()
当我们要为一个方阵M[m * m]
申请GPU的存储单元时
cudaMalloc((void **) &d_m,sizeof(int) * m * m);
d_m:指向存储在Device端数据的地址的指针(双重指针)。
sizeof(int) * m * m:存储在Device端空间的大小。
当我们要释放申请的GPU的存储单元时
cudaFree(d_m);
d_m:指向存储在Device端数据的地址的指针。
内存拷贝
cudaMemcpy(void dst, const void src, size_t count, cudaMemcpyKind kind)
- dst: 目的内存地址
- src: 源内存地址
- count: 以字节为单位要拷贝的数据大小
- kind: 数据传输方向
cudaMemcpyKind
- cudaMemcpyHostToDevice
- cudaMemcpyDeviceToHost
- cudaMemcpyDeviceToDevice
- cudaMemcpyHostToHost
当我们要将准备好的数据从CPU内存传输到GPU的存储单元时
cudaMemcpy(d_m, h_m, sizeof(int) * m * m, cudaMemcpyHostToDevice);
d_m:传输的目的地,GPU存储单元
h_m:数据的源地址,CPU存储单元
sizeof(int) * m * m:数据传输的大小
cudaMemcpyHostToDevice:数据传输的方向,CPU -> GPU
当我们要将准备好的数据从GPU存储单元传输到CPU的内存时
cudaMemcpy(h_c, d_c, sizeof(int) * m * k, cudaMemcpyDeviceToHost);
实验(矩阵乘法)
矩阵相乘计算样例:
CPU矩阵乘法
使用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.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;
}
}
}
可见整个计算过程需要完成三重循环,总共需要计算m * n * k
次tmp += h_m[i * n + h] * h_n[h * k + j];
才能完成计算。下面采用CUDA的并行编程模型来尝试去掉部分循环进行加速。
二维网格和线程块
为了解决矩阵乘法的问题,需要先搞定二维状态下线程的索引计算(当然也可以转换为一维)。
线程的索引排序方式如下图:
Global linear memory index: idx = iy * nx + ix
以44号线程为例
int index = (blockIdx.y * blockDim.y + threadIdx.y) * nx+ (blockIdx.x * blockDim.x + threadIdx.x);
Index = (2 * 2 + 1) * 8 + (1 * 4 + 0)
Index = 44
将整个grid对应到整个矩阵后,就如下图所示:
GPU矩阵乘法
grid中多线程计算Pd
- 每个线程计算Pd的一个元素
每个线程
- 读入矩阵Md的一行,读入矩阵Nd的一列
- 为每对Md和Nd元素执行一次乘法和加法
完整代码
#include
#include
#define BLOCK_SIZE 16
__global__ void gpu_matrix_mult(int *a,int *b, int *c, int m, int n, int k)
{
int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;
int sum = 0;
if( col < k && row < m)
{
for(int i = 0; i < n; i++)
{
sum += a[row * n + i] * b[i * k + col];
}
c[row * k + col] = sum;
}
}
void cpu_matrix_mult(int *h_a, int *h_b, 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.0;
for (int h = 0; h < n; ++h)
{
tmp += h_a[i * n + h] * h_b[h * k + j];
}
h_result[i * k + j] = tmp;
}
}
}
int main(int argc, char const *argv[])
{
int m=1000;
int n=2000;
int k=3000;
int *h_a, *h_b, *h_c, *h_cc;
cudaMallocHost((void **) &h_a, sizeof(int)*m*n);
cudaMallocHost((void **) &h_b, sizeof(int)*n*k);
cudaMallocHost((void **) &h_c, sizeof(int)*m*k);
cudaMallocHost((void **) &h_cc, sizeof(int)*m*k);
for (int i = 0; i < m; ++i) {
for (int j = 0; j < n; ++j) {
h_a[i * n + j] = rand() % 1024;
}
}
for (int i = 0; i < n; ++i) {
for (int j = 0; j < k; ++j) {
h_b[i * k + j] = rand() % 1024;
}
}
int *d_a, *d_b, *d_c;
cudaMalloc((void **) &d_a, sizeof(int)*m*n);
cudaMalloc((void **) &d_b, sizeof(int)*n*k);
cudaMalloc((void **) &d_c, sizeof(int)*m*k);
// copy matrix A and B from host to device memory
cudaMemcpy(d_a, h_a, sizeof(int)*m*n, cudaMemcpyHostToDevice);
cudaMemcpy(d_b, h_b, sizeof(int)*n*k, cudaMemcpyHostToDevice);
unsigned int grid_rows = (m + BLOCK_SIZE - 1) / BLOCK_SIZE;
unsigned int grid_cols = (k + BLOCK_SIZE - 1) / BLOCK_SIZE;
dim3 dimGrid(grid_cols, grid_rows);
dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);
gpu_matrix_mult<<>>(d_a, d_b, d_c, m, n, k);
cudaMemcpy(h_c, d_c, sizeof(int)*m*k, cudaMemcpyDeviceToHost);
cpu_matrix_mult(h_a, h_b, h_cc, m, n, k);
int ok = 1;
for (int i = 0; i < m; ++i)
{
for (int j = 0; j < k; ++j)
{
if(fabs(h_cc[i*k + j] - h_c[i*k + j])>(1.0e-10))
{
ok = 0;
}
}
}
if(ok)
{
printf("Pass!!!\n");
}
else
{
printf("Error!!!\n");
}
// free memory
cudaFree(d_a);
cudaFree(d_b);
cudaFree(d_c);
cudaFreeHost(h_a);
cudaFreeHost(h_b);
cudaFreeHost(h_c);
return 0;
}
代码解析
核函数
__global__ void gpu_matrix_mult(int *a,int *b, int *c, int m, int n, int k)
{
int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;
int sum = 0;
if( col < k && row < m)
{
for(int i = 0; i < n; i++)
{
sum += a[row * n + i] * b[i * k + col];
}
c[row * k + col] = sum;
}
}
先计算好元素在矩阵中的行和列,每一个线程处理了矩阵Md的一行和矩阵Nd的一列,分别相乘再相加。将计算好的最终值写入Pd矩阵的对应位置即可。
这里需要注意索引的计算方式,我们输入的矩阵实际上是一个一维的数组,那么需要将二维坐标线性化。
CPU计算函数
void cpu_matrix_mult(int *h_a, int *h_b, 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.0;
for (int h = 0; h < n; ++h)
{
tmp += h_a[i * n + h] * h_b[h * k + j];
}
h_result[i * k + j] = tmp;
}
}
}
将CPU计算结果作为标准来对GPU计算结果进行校验。
主函数
int *h_a, *h_b, *h_c, *h_cc;
cudaMallocHost((void **) &h_a, sizeof(int)*m*n);
cudaMallocHost((void **) &h_b, sizeof(int)*n*k);
cudaMallocHost((void **) &h_c, sizeof(int)*m*k);
cudaMallocHost((void **) &h_cc, sizeof(int)*m*k);
定义变量并在Host上分配数组空间。
for (int i = 0; i < m; ++i) {
for (int j = 0; j < n; ++j) {
h_a[i * n + j] = rand() % 1024;
}
}
for (int i = 0; i < n; ++i) {
for (int j = 0; j < k; ++j) {
h_b[i * k + j] = rand() % 1024;
}
}
在Host端初始化数组。
int *d_a, *d_b, *d_c;
cudaMalloc((void **) &d_a, sizeof(int)*m*n);
cudaMalloc((void **) &d_b, sizeof(int)*n*k);
cudaMalloc((void **) &d_c, sizeof(int)*m*k);
定义变量,并调用cudaMalloc
在GPU端为变量分配内存空间。
cudaMemcpy(d_a, h_a, sizeof(int)*m*n, cudaMemcpyHostToDevice);
cudaMemcpy(d_b, h_b, sizeof(int)*n*k, cudaMemcpyHostToDevice);
将在CPU端初始化好的数组拷贝到GPU,便于在GPU并行处理数据。
dim3 dimGrid(grid_cols, grid_rows);
dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);
gpu_matrix_mult<<>>(d_a, d_b, d_c, m, n, k);
cudaMemcpy(h_c, d_c, sizeof(int)*m*k, cudaMemcpyDeviceToHost);
完成执行设置并调用核函数,在GPU上开始程序的并行部分。
计算完成后结果存储在d_c
所指向的GPU内存中,需要调用cudaMemcpy
将数据从Device拷贝到Host中。
cpu_matrix_mult(h_a, h_b, h_cc, m, n, k);
int ok = 1;
for (int i = 0; i < m; ++i)
{
for (int j = 0; j < k; ++j)
{
if(fabs(h_cc[i*k + j] - h_c[i*k + j])>(1.0e-10))
{
ok = 0;
}
}
}
if(ok)
{
printf("Pass!!!\n");
}
else
{
printf("Error!!!\n");
}
在CPU中完成标准答案的计算并与GPU计算结果进行对比校验。
cudaFree(d_a);
cudaFree(d_b);
cudaFree(d_c);
cudaFreeHost(h_a);
cudaFreeHost(h_b);
cudaFreeHost(h_c);
调用cudaFree
相关函数释放刚刚在Host和Device中分配的内存,随即程序结束。