本系列代码托管于:https://github.com/chintsan-code/cuda_by_example
本篇使用的项目为:add_loop_blocks、add_loop_long_blocks、add_loop_very_long_blocks
在之前的篇章(CUDA实战笔记(4)——矢量求和运算)中,在基于GPU实现矢量求和部分,我们通过add<<N, 1>>
告诉CUDA Runtime启动N个核函数副本,我们将这些副本称为线程块(Block)。这就是尖括号中第一个参数的作用。第二个参数表示CUDA Runtime在每个Block中创建的线程(Thread)数量,在那个例子例子中,每个Block中只包含一个Thread。因此add<N, 1>
总共让CUDA Runtime启动了N*1 = N个并行线程。
一. 使用线程实现GPU上的矢量求和
下面将使用Thread代替Block实现矢量求和。在这一问题上这种改变几乎没有任何优势,但是Block中的并行Thread能够实现并行Block无法完成的工作。
首先修改核函数add
的调用部分,之前是启动N个Block,每个Block包含1个Thread。现在要改为启动1个Block,它包含N个Thread:
add<<<1, N>>>(dev_a, dev_b, dev_c);
接着要修改核函数add
,之前是通过Block的索引blockIdx
实现数输入输出的数据进行索引的,现在只启动1个Block,自然是要对着部分进行修改:
__global__ void add(int* a, int* b, int* c)
{
int tid = threadIdx.x; // 只有一个Block,可以通过线程索引来对数据进行索引
if (tid < N)
{
c[tid] = a[tid] + b[tid];
}
}
这样就完成了修改,完整代码如下:
// add_loop_blocks
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <stdio.h>
#include "../../common/book.h"
#define N 10
__global__ void add(int* a, int* b, int* c) {
int tid = threadIdx.x; // 只有一个Block,通过线程索引来对数据进行索引
if (tid < N) {
c[tid] = a[tid] + b[tid];
}
}
int main() {
int a[N], b[N], c[N];
int* dev_a, * dev_b, * dev_c;
HANDLE_ERROR(cudaMalloc((void**)&dev_a, N * sizeof(int)));
HANDLE_ERROR(cudaMalloc((void**)&dev_b, N * sizeof(int)));
HANDLE_ERROR(cudaMalloc((void**)&dev_c, N * sizeof(int)));
for (int i = 0; i < N; i++)
{
a[i] = -i;
b[i] = i * i;
}
HANDLE_ERROR(cudaMemcpy(dev_a, a, N * sizeof(int), cudaMemcpyHostToDevice));
HANDLE_ERROR(cudaMemcpy(dev_b, b, N * sizeof(int), cudaMemcpyHostToDevice));
add<<<1, N>>>(dev_a, dev_b, dev_c); // 启动1个Block,它包含N个Thread
HANDLE_ERROR(cudaMemcpy(c, dev_c, sizeof(int) * N, cudaMemcpyDeviceToHost));
for (int i = 0; i < N; i++)
{
printf("%d + %d = %d\n", a[i], b[i], c[i]);
}
HANDLE_ERROR(cudaFree(dev_a));
HANDLE_ERROR(cudaFree(dev_b));
HANDLE_ERROR(cudaFree(dev_c));
return 0;
}
下面会对这种只使用Thread的方法中存在的局限进行说明:
二. 在GPU上对更长的矢量求和
在《CUDA实战笔记(4)——矢量求和运算》的末尾,有提到由于硬件限制,Block的数量是有限的,同样的,每一个Block中的Thread数量也是有限的。可以通过cudaDeviceProp.maxThreadsPerBlock
获取每一个Block中的最大Thread数量,我这里是1024。
现在使用多个Block,每个Block包含多个Thread,计算索引看上去非常类似于将二维索引空间转换为一维空间:
int tid = blockIdx.x * blockDim.x + threadIdx.x;
这里出现了一个全新的内置变量blockDim
。对于所有Block来说,这个变量是一个常数,它等于Block中每一维的Thread的数量。由于本次修改只用到一维的Block,因此只用到blockDim.x
,类似之前使用过的gridDim
。本次问题的线程块和线程的的空间组织形式如下:
如果用Block表示行,用Thread表示列,那么就可以计算得到唯一的索引。例如要计算上图中绿色背景的Thread在在线性空间中的唯一索引,那么就是:
int tid = blockIdx.x * blockDim.x + threadIdx.x; // 2 * 4 + 2 = 10
既然修改为多个Block和多个Thread的模型,那么调用核函数的部分也需要一起进行修改,仍然是启动N个并行线程,但是为了不超过Block中Thread的最大限制1024,将每一个Block中的Thread数量设置为128,然后启动N/128个Block,这样就启动了N个Thread同时并行运行。
但是这样会引入一个新的问题,如果N不是128的倍数,例如N为127,那么就会启动0个Block,这就无法计算,实际上,如果N不是128的倍数,启动的Block都会少于预期数量,因此我们希望N/128的结果可以向上取整:
add<<<(127 + N) / 128, 128>>>(dev_a, dev_b, dev_c); // (127 + N) / 128 : N/128向上取整
上面的代码保证了当N不为128的整数倍时,启动比预期更多的Block,既然有可能启动多余的Block和Thread,那么久需要有一种机制,可以保证这些多余的Block和Thread不参与运算,其实之前我们一直在使用这种机制:
int tid = blockIdx.x * blockDim.x + threadIdx.x;
if (tid < N) {
c[tid] = a[tid] + b[tid];
}
当索引tid越过数组时,核函数会自动停止执行计算,更重要的是,核函数不会越过数组边界对内存进行读写。
完整代码如下:
// add_loop_long_blocks
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <stdio.h>
#include "../../common/book.h"
#define N (33*1024) // 这么大的N,如果只用一个Block是做不了的
__global__ void add(int* a, int* b, int* c) {
int tid = blockIdx.x * blockDim.x + threadIdx.x;
if (tid < N) {
c[tid] = a[tid] + b[tid];
}
}
// 当N很大时要修改项目属性->配置属性->链接器->系统->堆栈保留大小(单位: byte)
// 或者分配在全局区或堆区
// int a[N], b[N], c[N]; // 分配在全局区
int main() {
//int a[N], b[N], c[N]; // 分配在栈区
int* a, * b, * c;
a = (int*)malloc(N * sizeof(int)); // 分配在堆区,最后记得要释放
b = (int*)malloc(N * sizeof(int));
c = (int*)malloc(N * sizeof(int));
int* dev_a, * dev_b, * dev_c;
HANDLE_ERROR(cudaMalloc((void**)&dev_a, N * sizeof(int)));
HANDLE_ERROR(cudaMalloc((void**)&dev_b, N * sizeof(int)));
HANDLE_ERROR(cudaMalloc((void**)&dev_c, N * sizeof(int)));
for (int i = 0; i < N; i++)
{
a[i] = -i;
b[i] = i * i;
}
HANDLE_ERROR(cudaMemcpy(dev_a, a, N * sizeof(int), cudaMemcpyHostToDevice));
HANDLE_ERROR(cudaMemcpy(dev_b, b, N * sizeof(int), cudaMemcpyHostToDevice));
// (127 + N) / 128 : N/128向上取整
// 每个Block有128个Thread
// 注意( N + 127 ) / 128不能超过maxGridSize的限制
add<<<(127 + N) / 128, 128>>>(dev_a, dev_b, dev_c);
HANDLE_ERROR(cudaMemcpy(c, dev_c, sizeof(int) * N, cudaMemcpyDeviceToHost));
for (int i = 0; i < N; i++)
{
printf("%d + %d = %d\n", a[i], b[i], c[i]);
}
HANDLE_ERROR(cudaFree(dev_a));
HANDLE_ERROR(cudaFree(dev_b));
HANDLE_ERROR(cudaFree(dev_c));
free(a); // 释放之前分配的堆区内存
free(b);
free(c);
return 0;
}
这样修改之后基本可以满足相当大的N长度的矢量相加了,但是受限于硬件,Grid中包含的Block也是有限的,可以通过cudaDeviceProp.maxGridSize
获取Grid每一维度的Block上限,我这里是(2147483647, 65535, 65535),因此如果上面的N超过了(2^31 – 1)*128,那么核函数就会调用失败(当然在调用核函数之前,由于显存不足就会报错:out of memory了)。
假设我们的Grid最大可以包含65535个Block(以前的GPU是这样的),且不会出现显存不足的问题。那么当如果上面的N超过了 65535 *128,那么核函数就会调用失败。下面来解决这个问题。
三. 在GPU上对任意长度的矢量求和
首先对核函数进行一些修改:
__global__ void add(int* a, int* b, int* c) {
int tid = blockIdx.x * blockDim.x + threadIdx.x;
while (tid < N) {
c[tid] = a[tid] + b[tid];
tid += gridDim.x * blockDim.x;
}
}
这类似于最初CPU版本的实现:
void add(int* a, int* b, int* c)
{
int tid = 0; // 这是第0个CPU,因此索引从0开始
while (tid < N)
{
c[tid] = a[tid] + b[tid];
tid += 1; // 由于我们只有一个CPU,因此每次递增1
}
}
在这里,使用了一个while循环对数据进行迭代,在CPU版本中,每次递增的数量是CPU的数量,在GPU中的实现也是同样的道理:由于硬件的限制,有时候没办法让每个数据都被一个Thread独立处理,因此可以让一个Thread处理多个数据。
我们希望每个并行Thread从不同的索引开始,因此每个线程的索引按照以下公式获取:
int tid = blockIdx.x * blockDim.x + threadIdx.x;
当该线程处理完一个数据之后,加上一个步长,处理下一个数据,直到数据的索引超过矢量长度N。关键点就在这个步长应该怎么得到。可以这么理解:当启动了10个Thread,要处理100个数据,而1个Thread只能同时处理一个数据。那么有:
- Thread0处理:第0、10、20、30、40、50、60、70、80、90个数据
- Thread1处理:第1、11、21、31、41、51、61、71、81、91个数据
- …
- Thread9处理:第9、19、29、39、49、59、69、79、89、99个数据
可以看到,步长其实就等于启动的Thread数量,而
$$启动的Thread数量 = 启动的Block数量 * 每个Block包含的线程数量$$
于是有
tid += gridDim.x * blockDim.x;
其中gridDim.x
就是在这个矢量求和问题中启动的Block数量,blockDim.x
就是每个Block中包含Thread的数量。
核函数就修改完成了,接下来修改它的调用部分。前面说过,如果启动的Block大于硬件设置,就会调用失败,因此要减少Block的数量,这里将Block的数量固定为某一个较小的值。例如128:
add<<<128, 128 >>>(dev_a, dev_b, dev_c);
这样我们就启动了128个Block,每个Block包含128个Thread。
在之前,总共启动N个Thread,每个Thread只需要处理1个数据。现在是启动128*128个Thread,每个Thread分别要处理N/16384个数据。
当前,也可以调节Block和Thread到合适的大小。
下面是完整代码。
// add_loop_very_long_blocks
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <stdio.h>
#include "../../common/book.h"
#define N (65536*1024)
__global__ void add(int* a, int* b, int* c) {
int tid = blockIdx.x * blockDim.x + threadIdx.x;
while (tid < N) {
c[tid] = a[tid] + b[tid];
tid += gridDim.x * blockDim.x; // 步长:启动的Thread的数量 = 启动的Block数量(gridDim.x) * 每个Block包含的Thread的数量(blockDim.x)
// 让每个Thread处理 N / (gridDim.x * blockDim.x)个数据
}
}
int main() {
int* a, * b, * c;
a = (int*)malloc(N * sizeof(int));
b = (int*)malloc(N * sizeof(int));
c = (int*)malloc(N * sizeof(int));
int* dev_a, * dev_b, * dev_c;
HANDLE_ERROR(cudaMalloc((void**)&dev_a, N * sizeof(int)));
HANDLE_ERROR(cudaMalloc((void**)&dev_b, N * sizeof(int)));
HANDLE_ERROR(cudaMalloc((void**)&dev_c, N * sizeof(int)));
for (int i = 0; i < N; i++)
{
a[i] = -i;
b[i] = i * i;
}
HANDLE_ERROR(cudaMemcpy(dev_a, a, N * sizeof(int), cudaMemcpyHostToDevice));
HANDLE_ERROR(cudaMemcpy(dev_b, b, N * sizeof(int), cudaMemcpyHostToDevice));
// 启动128个Block,每个Block含有128个线程
// 让每个Thread处理 N / 16384个数据
add<<<128, 128>>>(dev_a, dev_b, dev_c);
HANDLE_ERROR(cudaMemcpy(c, dev_c, sizeof(int) * N, cudaMemcpyDeviceToHost));
for (int i = 0; i < N; i++)
{
printf("%d + %d = %d\n", a[i], b[i], c[i]);
}
HANDLE_ERROR(cudaFree(dev_a));
HANDLE_ERROR(cudaFree(dev_b));
HANDLE_ERROR(cudaFree(dev_c));
free(a);
free(b);
free(c);
return 0;
}
参考:
- 《GPU高性能编程 CUDA实战》
评论 (0)