CUDA实战笔记(13)——原子性

本系列代码托管于：https://github.com/chintsan-code/cuda_by_example

本篇使用的项目为：hist_cpu、hist_gpu_gmem_atomics、hist_gpu_shmem_atomics

首先来看一个简单的操作：

x++;

这个简单的语句包含3个操作：

1. 读取x中的值；
2. 将步骤1中读取到的值增加1；
3. 将步骤2得到的结果写回变量x。

这个过程也被称为读取-修改-写入(Read-Modify-Write)操作，其中步骤2的递增操作也可以换成其他修改x值额操作。

假如现在有2个线程都需要对x的值进行上述修改，这两个线程用A和B表示，假设x的初始值为7，那么理想的情况是下面这种，此时会得到正确的结果：

1. 线程A读取x中的值(x=7)
2. 线程A将读取到的值增加1(7+1=8)
3. 线程A将结果写回x(x=8)
4. 线程A读取x中的值(x=8)
5. 线程A将读取到的值增加1(8+1=9)
6. 线程A将结果写回x(x=9)

但是，如果线程调度方式有差异，也可能出现下面这种异常情况，得到错误的结果：

1. 线程A读取x中的值(x=7)
2. 线程B读取x中的值(x=7)
3. 线程A将读取到的值增加1(7+1=8)
4. 线程B将读取到的值增加1(7+1=8)
5. 线程A将结果写回x(x=8)
6. 线程B将结果写回x(x=8)

因此，需要有某种方式，可以一次性地执行读取-修改-写入(Read-Modify-Write)这三个操作，除非已经完成了这三个操作，否则其他线程都不能读取或者写入x的值。(也就是加锁)

由于上面这种方式的执行过程不能分解为更小部分，因此把满足这种条件限制的操作称为原子操作。

下面通过一个计算直方图的例子，演示CUDA C的原子操作。

一. 在CPU上计算直方图

为了进行性能测试，需要有一个基准性能，这里编写了一个CPU程序：

// hist_cpu

#include "../../common/book.h"
#include "time.h"

#define SIZE (100*1024*1024)

int main() {
    // 随机生成100MB的随机数据
    unsigned char* buffer = (unsigned char*)big_random_block(SIZE);
    // 每个字节的取值范围为0x00-0xFF,因此使用大小为256的数组来存储相应的值在buffer中出现的次数,
    // 用于计算直方图
    unsigned int histo[256] = { 0 };
    // 计算时间
    clock_t start, stop;
    start = clock();
    for (int i = 0; i < SIZE; i++) {
        histo[buffer[i]]++;
    }
    stop = clock();
    float elapsedTime = (float)(stop - start) / (float)CLOCKS_PER_SEC * 1000.0f;
    printf("Time to generate:  %3.1f ms\n", elapsedTime);
    // 验证直方图的所有元素加起来是否等于正确的值(应该等于SIZE)
    long histoCount = 0;
    for (int i = 0; i < 256; i++) {
        histoCount += histo[i];
    }
    printf("Histogram Sum:  %ld\n", histoCount);

    // 释放内存
    free(buffer);

    return 0;
}

(注意：这里没有使用多线程)

计算一下处理时间：46.0 ms

二. 在GPU上计算直方图

现在，使用CUDA计算直方图，首先完成除了核函数部分的代码：

int main() {
    // 随机生成100MB的随机数据
    unsigned char* buffer = (unsigned char*)big_random_block(SIZE);

    // 初始化计时事件
    cudaEvent_t start, end;
    HANDLE_ERROR(cudaEventCreate(&start));
    HANDLE_ERROR(cudaEventCreate(&end));
    HANDLE_ERROR(cudaEventRecord(start, 0));

    // 在GPU上为文件的数据分配内存
    unsigned char* dev_buffer;
    unsigned int* dev_histo;
    HANDLE_ERROR(cudaMalloc((void**)&dev_buffer, SIZE * sizeof(char)));
    HANDLE_ERROR(cudaMalloc((void**)&dev_histo, 256 * sizeof(int)));
    HANDLE_ERROR(cudaMemcpy(dev_buffer, buffer, SIZE * sizeof(char), cudaMemcpyHostToDevice));
    HANDLE_ERROR(cudaMemset(dev_histo, 0, 256 * sizeof(int)));

    cudaDeviceProp prop;
    HANDLE_ERROR(cudaGetDeviceProperties(&prop, 0));
    int blocks = prop.multiProcessorCount * 2; // 将Block的数量设置为GPU中处理器数量的2倍
    histo_kernel<<<blocks, 256>>>(dev_buffer, SIZE, dev_histo);

    unsigned int histo[256];
    HANDLE_ERROR(cudaMemcpy(histo, dev_histo, 256 * sizeof(int), cudaMemcpyDeviceToHost));

    // 得到停止时间并显示计时结果
    HANDLE_ERROR(cudaEventRecord(end, 0));
    HANDLE_ERROR(cudaEventSynchronize(end));
    float elapsedTime;
    HANDLE_ERROR(cudaEventElapsedTime(&elapsedTime, start, end));
    printf("Time to generate:  %3.1f ms\n", elapsedTime);

    // 验证直方图的所有元素加起来是否等于正确的值(应该等于SIZE)
    long histoCount = 0;
    for (int i = 0; i < 256; i++) {
        histoCount += histo[i];
    }
    printf("Histogram Sum:  %ld\n", histoCount);

    // 验证与CPU得到的是相同的计数值
    for (int i = 0; i < SIZE; i++) {
        histo[buffer[i]]--;
    }
    for (int i = 0; i < 256; i++) {
        if (histo[i] != 0) {
            printf("Failure at %d!  Off by %d\n", i, histo[i]);
        }
    }

    // 释放事件、内存
    HANDLE_ERROR(cudaEventDestroy(start));
    HANDLE_ERROR(cudaEventDestroy(end));
    HANDLE_ERROR(cudaFree(dev_buffer));
    HANDLE_ERROR(cudaFree(dev_histo));
    free(buffer);

    return 0;
}

由于本次问题的直方图包含256个元素，因此把每个Block中的Thread数量设置为256。

这里为了验证在GPU上计算的直方图结果是正确的，在程序的后半部分，使用CPU来计算直方图作为验证手段。

接下来完成核函数部分。

三. 使用全局内存原子操作的直方图核函数

首先要将获取Grid中的线程ID：

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

由于在上面设计线程模型时，每个Thread都要不止计算一个数据，因此在一个Thread计算完之后，需要接着计算下一个数据：

int stride = gridDim.x * blockDim.x;
while (tid < size) {
    // TODO: 计算直方图
    tid += stride;
}

这里的步长就是Grid中的线程数量，假设Grid中只有1个Block，Block中有256个线程(即步长stride为1*256=256)，那么有：

• Thread0：计算buffer[0]、buffer[256]、buffer[512]…
• Thread1：计算buffer[1]、buffer[257]、buffer[513]…
• ……
• Thread255：计算buffer[255]、buffer[511]、buffer[767]…

接下来需要一个原子操作，因为虽然在while中，同一个Thread只有当上一次数据计算完成之后才会执行下一个，属于串行计算。但是Block中有256个Thread，因此它们可能会同时访问histo[]中的同一个位置，因此需要一个原子操作：

atomicAdd(&(histo[buffer[tid]]), 1);

这里表示在CUDA C中使用原子操作的方式。函数调用atomicAdd(address, val)将生成一个原子的序列操作，这个操作序列包括读取地址address处的值，将val增加到这个值上，以及将结果保存回地址address。底层硬件将确保当执行这些操作时，其他任何线程都不会读取或写入地址address上的值，这样就能确保得到预计的结果。

假如这里没有使用原子操作，而是：

histo[buffer[tid]]++;

那么，如果buffer[0]和buffer[1]的数值相同，如果线程调度出现问题，就可能出现错误的结果：

1. Thread0读取buffer[0]，这里的值假设为a(0≤a≤255)
2. Thread1读取buffer[1] ，这里的值同样为a
3. Thread0读取histo[a]，由于之前没有计算，因此值为0
4. Thread1读取histo[a]，由于之前没有计算，因此值为0
5. Thread0计算直方图：histo[a]++，得到结果1
6. Thread1计算直方图：histo[a]++，得到结果1

这与正确结果histo[a]等于2不一致。因此在GPU上计算直方图，要使用原子操作。

完整的核函数代码如下：

__global__ void histo_kernel(unsigned char* buffer,
                             long size, 
                             unsigned int* histo) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    int stride = gridDim.x * blockDim.x;
    while (tid < size) {
        atomicAdd(&(histo[buffer[tid]]), 1);
        tid += stride;
    }
}

这样就完成了GPU上的直方图，计算处理时间：74.1 ms。

这里甚至比不使用多线程的CPU程序耗时还要长。这是因为全局内存上的原子操作降低了程序性能：当数千个Thread尝试同时访问少量的内存位置时(这里就是访问histo[]，只有256个位置)，将发生大量的竞争，这将降低性能。

下面对程序算法进行改进。

四. 使用共享内存原子操作和全局内存原子操作的直方图核函数

上面分析了性能下降的原因是当数千个Thread尝试同时访问少量的内存位置时(这里就是访问histo[]，只有256个位置)，将发生大量的竞争。因此改进的思路就是降低同时访问内存位置的线程数量。

具体方法是：将这一计算分为两个阶段来进行：第一阶段，每一个Block都创建一个共享缓冲区，用来计算这个Block中对应buffer中的数据的直方图；第二阶段，当所有Block都计算完成之后，再把它们相加起来。

首先分配一个共享内存缓冲区并进行初始化，用来保存每个Block的临时直方图：

__shared__ unsigned int temp[256];
temp[threadIdx.x] = 0;
__syncthreads();

注意这里要加个个同步线程的操作，只有这样，才可以包装同一个Block中temp[]的所有256个元素都被初始化为0了。

接着和全局原子操作类似，只是这里将histo[]改成了temp[]：

int tid = blockIdx.x * blockDim.x + threadIdx.x;
int stride = gridDim.x * blockDim.x;
while (tid < size) {
    atomicAdd(&(temp[buffer[tid]]), 1);
    tid += stride;
}

假设这个Grid中有4个Block，每个Block中有256个Thread(即步长stride为4*256=1024)，那么会创建4个共享缓冲区，为了方便说明，用temp0[]，temp1[]，temp2[]，temp3[]来表示。那么有：

在Block0中：

• Thread0：计算buffer[0]、buffer[1024]、buffer[2048]…将结果暂时保存到Block0中的共享缓冲区temp0[]
• Thread1：计算buffer[1]、buffer[1025]、buffer[2049]… 将结果暂时保存到Block0中的共享缓冲区temp0[]
• ……
• Thread255：计算buffer[255]、buffer[1279]、buffer[2303]… 将结果暂时保存到Block0中的共享缓冲区temp0[]

在Block1中：

• Thread256：计算buffer[256]、buffer[1280]、buffer[2304]…将结果暂时保存到Block1中的共享缓冲区temp1[]
• Thread257：计算buffer[257]、buffer[1281]、buffer[2305]… 将结果暂时保存到Block1中的共享缓冲区temp1[]
• ……
• Thread511：计算buffer[511]、buffer[1535]、buffer[2559]… 将结果暂时保存到Block1中的共享缓冲区temp1[]

在Block2中：

• Thread512：计算buffer[512]、buffer[1536]、buffer[2560]…将结果暂时保存到Block2中的共享缓冲区temp2[]
• Thread513：计算buffer[257]、buffer[1537]、buffer[2561]… 将结果暂时保存到Block2中的共享缓冲区temp2[]
• ……
• Thread767：计算buffer[767]、buffer[1791]、buffer[2815]… 将结果暂时保存到Block2中的共享缓冲区temp2[]

在Block3中：

• Thread768：计算buffer[768]、buffer[1792]、buffer[2816]…将结果暂时保存到Block3中的共享缓冲区temp3[]
• Thread769：计算buffer[769]、buffer[1793]、buffer[2817]… 将结果暂时保存到Block3中的共享缓冲区temp3[]
• ……
• Thread1023：计算buffer[1023]、buffer[2047]、buffer[3071]… 将结果暂时保存到Block3中的共享缓冲区temp3[]

之后在汇总所有共享缓冲区，即可得到最后的直方图：

__syncthreads();
atomicAdd(&histo[threadIdx.x], temp[threadIdx.x]);

注意这里也要同步线程，因为要等待所有的Block都完成计算。

完整的核函数代码如下：

__global__ void histo_kernel(unsigned char* buffer,
                             long size,
                             unsigned int* histo) {
    // 分配一个共享内存缓冲区并进行初始化，用来保存每个线程块的临时直方图
    __shared__ unsigned int temp[256];
    temp[threadIdx.x] = 0;
    __syncthreads();  // 确保每个线程的写入操作(temp初始化)在线程继续前进之前完成

    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    int stride = gridDim.x * blockDim.x;
    while (tid < size) {
        atomicAdd(&(temp[buffer[tid]]), 1);
        tid += stride;
    }
    __syncthreads();
    atomicAdd(&histo[threadIdx.x], temp[threadIdx.x]);
}

计算处理时间：20.4 ms

---

参考：

• 《GPU高性能编程 CUDA实战》