CUDA实战笔记(6)——使用线程实现GPU上的矢量求和萤火

本系列代码托管于：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实战》

CUDA实战笔记(6)——使用线程实现GPU上的矢量求和

一. 使用线程实现GPU上的矢量求和

二. 在GPU上对更长的矢量求和

三. 在GPU上对任意长度的矢量求和

评论 (0)

专题

CUDA编程入门¹⁴

Linux系统编程¹²

OpenCV4入门³⁷

机器学习入门¹⁸

相关文章

CUDA实战笔记(6)——使用线程实现GPU上的矢量求和

一. 使用线程实现GPU上的矢量求和

二. 在GPU上对更长的矢量求和

三. 在GPU上对任意长度的矢量求和

评论 (0)

猜你喜欢

专题

CUDA编程入门14

Linux系统编程12

OpenCV4入门37

机器学习入门18

相关文章

CUDA编程入门¹⁴

Linux系统编程¹²

OpenCV4入门³⁷

机器学习入门¹⁸