CUDA By Example(七)——原子性

最新推荐文章于 2025-05-06 14:41:11 发布

hhhcbw

最新推荐文章于 2025-05-06 14:41:11 发布

阅读量870

点赞数

分类专栏： GPU高性能编程CUDA实战文章标签： c++ CUDA 多线程并行计算 GPU高性能编程

本文链接：https://blog.youkuaiyun.com/weixin_44491423/article/details/128906714

版权

GPU高性能编程CUDA实战专栏收录该内容

8 篇文章

订阅专栏

文章目录

原子操作简介
计算直方图
更多相关内容

原子操作简介

在编写传统的单线程应用程序时，程序员通常不需要使用原子操作。即使需要使用原子操作也无需担心。下面将详细解释原子操作是什么，以及为什么在多线程程序中需要使用它们。为了解释原子操作，首先来分析 C 或者 C++ 的基础知识之一，即递增运算符：

x++;

这是在标准 C 中的一个表达式，在执行完这条语句后，x 的值应该比执行这条语句之前的值大 1。这条语句中包含如下三个步骤：

读取 x 中的值
将步骤 1 中读到的值增加 1
将递增后的结果写回到 x

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

现在考虑一种情况：有两个线程都需要对 x 的值进行递增。将这两个线程分别称为 A 和 B。A 和 B 在递增 x 的值时都需要执行上面三个操作。假设 x 的初始值为 7。在理想情况下，我们希望线程 A 和 B 执行下表中的步骤

步骤	示例
1) 线程 A 读取 `x` 中的值	A 从 `x` 中读到 7
2) 线程 A 将读到的值增加 1	A 计算得到 8
3) 线程 A 将结果写回到 `x`	`x` <- 8
4) 线程 B 读取 `x` 中的值	B 从 `x` 中读到 8
5) 线程 B 将读到的值增加 1	B 计算得到 9
6) 线程 B 将结果写回到 `x`	`x` <- 9

由于 x 的起始值为 7，并且由两个线程进行递增，因此在递增运算完成后，x 的值变为 9。根据前面的操作顺序，这确实是文明得到的结果。遗憾的是，除了这个操作顺序外，还存在其他一些操作顺序可能导致错误的结果。例如，下表中的顺序，其中线程 A 和 B 的操作彼此交叉进行。

步骤	示例
1) 线程 A 读取 `x` 中的值	A 从 `x` 中读到 7
2) 线程 B 读取 `x` 中的值	B 从 `x` 中读到 7
3) 线程 A 将读到的值增加 1	A 计算得到 8
4) 线程 B 将读到的值增加 1	B 计算得到 8
5) 线程 A 将结果写回到 `x`	`x` <- 8
6) 线程 B 将结果写回到 `x`	`x` <- 8

因此，如果线程的调度方式不正确，那么最终将得到错误的结果。除了上面两种执行顺序外，这 6 个步骤还有许多其他的排序方式，其中有些方式能产生正确的结果，而其他的方式则不能。当把程序从单线程改写为多线程时，如果多个线程需要对共享值进行读取或者写入时，那么很可能会遇到不可预测的结果。

在前面的示例中，我们需要通过某种方式一次性地执行完读取-修改-写入这三个操作，并且在执行过程中不会被其他线程中断。具体来说，除非已经完成了这三个操作，否则其他的线程都不能读取或者写入 x 的值。由于这些操作的执行过程不能分解为更小的部分，因此我们将满足这种条件限制的操作称为原子操作。CUDA C 支持多种原子操作，当有数千个线程在内存访问上发生竞争时，这些操作能够确保在内存上实现安全的操作。

计算直方图

给定一个包含一组元素的数据集，直方图表示每个元素出现的频率。例如，如果计算词组 “Programming with CUDA C” 中字符频率的直方图，那么将得到下表结果。

A	C	D	G	H	I	M	N	O	P	R	T	U	W
2	2	1	2	1	2	2	1	1	1	2	1	1	1

虽然直方图的定义很简单，但却在计算机科学领域得到了非常广的应用。在各种算法中都用到直方图，包括图像处理、数据压缩、计算机视觉、机器学习、音频编码等等。下面将把直方图运算作为代码示例的算法。

在 CPU 上计算直方图

下面首先给出如何在CPU上计算直方图。这个示例同时也说明了在单线程应用中计算直方图是非常简单的。这个应用程序将处理一个大型的数据流。在实际程序中，这个数据可以是像素的颜色值，或者音频采样数据，但在我们的示例程序中，这个数据是随机生成的字节流。我们可以通过工具函数 big_random_block() 来生成这个随机的字节流。在应用程序中将生成 100MB 的随机数据。

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

#include <stdio.h>
#include <iostream>
#include <windows.h>

#define SIZE (100*1024*1024)

int main(void) {
	unsigned char* buffer = (unsigned char*)big_random_block(SIZE);

由于每个随机字节 (8比特) 都有 256 个不同的可能取值（从0x00到0xFF），因此在直方图中需要包含 256 个元素，每个元素记录相应的值在数据流中出现次数。

下面创建一个包含 256 个元素的数组，并将所有元素的值初始化为 0。

unsigned int histo[256];
for (int i = 0; i < 256; i++)
	histo[i] = 0;

在创建了直方图并将数组元素初始化为 0 后，接下来需要计算每个值在 buffer[] 数据中的出现频率。

算法的思想是，每当在数组 buffer[] 中出现某个值 z 时，就递增直方图数组中索引为 z 的元素。这样就能计算出值 z 的出现次数。

如果当前看到的值为 buffer[i]，那么将递增数组中索引等于 buffer[i] 的元素。由于元素 buffer[i] 位于 histo[buffer[i]]，我们只需一行代码就可以递增相应的计数器。

histo[buffer[i]]++;

我们在简单的 for() 循环中对 buffer[] 每个元素执行这个操作。

for (int i = 0; i < SIZE; i++)
	histo[buffer[i]]++;

此时，我们已经计算完了输入数据的直方图。在实际的应用程序中，这个直方图可能作为下一个计算步骤的输入数据。但在这里的简单示例中，这就是要执行的所有工作，因此接下来将验证直方图的所有元素相加起来是否等于正确的值，然后结束程序。

long histoCount = 0;
for (int i = 0; i < 256; i++) {
	histoCount += histo[i];
}
std::cout << "Histogram Sum: " << histoCount << std::endl;

思考一下就会发现，无论输入数组的值是什么，这个和值总是相同的。每个元素都将统计相应数值的出现次数，因此所有这些元素值的总和就应该等于数组中元素的总数量。在示例中，这个值就等于 SIZE。

在执行完运算后需要释放内存并返回。

	free(buffer);
	return 0;
}

完整代码

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

#include <stdio.h>
#include <iostream>
#include <windows.h>

#define SIZE (100*1024*1024)

int main(void) {
	unsigned char* buffer = (unsigned char*)big_random_block(SIZE);
	DWORD start, end;
	start = GetTickCount();
	unsigned int histo[256];
	for (int i = 0; i < 256; i++)
		histo[i] = 0;
	for (int i = 0; i < SIZE; i++)
		histo[buffer[i]]++;
	end = GetTickCount();
	std::cout << "Time to generate: " << end - start << " ms" << std::endl;
	long histoCount = 0;
	for (int i = 0; i < 256; i++) {
		histoCount += histo[i];
	}
	std::cout << "Histogram Sum: " << histoCount << std::endl;
	free(buffer);
	return 0;
}

运行结果

在GPU上计算直方图

我们把这个直方图计算示例改在 GPU 上运行。如果输入的数组足够大，那么通过由多个线程来处理缓冲区的不同部分，将节约大量的计算时间。其中，由不同的线程来读取不同部分的输入数据是非常容易的。

在计算输入数组的直方图时存在一个问题，即多个线程可能同时对输出直方图的同一个元素进行递增。在这种情况下，我们需要通过原子的递增操作来避免之前提到线程执行顺序不同带来的问题。

main() 函数的开头与基于 CPU 的版本完全一样

int main( void ) {
    unsigned char* buffer = (unsigned char*)big_random_block(SIZE);

由于要测量代码的执行性能，因此要初始化计时事件。

cudaEvent_t start, stop;
HANDLE_ERROR(cudaEventCreate(&start));
HANDLE_ERROR(cudaEventCreate(&stop));
HANDLE_ERROR(cudaEventRecord(start, 0));

在设置好输入数据和事件后，我们需要在 GPU 上为随机输入数据和输出直方图分配内存空间。在分配了输入缓冲区后，我们将 big_random_block() 生成的数组复制到 GPU 上。同样，在分配了直方图后，像 CPU 版本中那样将其初始化为 0。

// 在GPU上为文件的数据分配内存
unsigned char* dev_buffer;
unsigned int* dev_histo;

HANDLE_ERROR(cudaMalloc((void**)&dev_buffer, SIZE));
HANDLE_ERROR(cudaMemcpy(dev_buffer, buffer, SIZE, cudaMemcpyHostToDevice));
HANDLE_ERROR(cudaMalloc((void**)&dev_histo, 256 * sizeof(int)));
HANDLE_ERROR(cudaMemset(dev_histo, 0, 256 * sizeof(int)));

cudaMemset() 这个函数的原型与标准 C 函数 memset() 的原型是相似的，并且这两个函数的行为也基本相同。二者的差异在于，cudaMemset() 将返回一个错误码，而 C 库函数 memset() 则不是。这个错误码将告诉调用者在设置 GPU 内存时发生的错误。除了返回错误码外，还有一个不同之处就是，cudaMemset() 是在 GPU 内存上执行，而 memset() 是在主机内存上运行。

在初始化输入缓冲区和输出缓冲区后，就做好了计算直方图的准备。你马上就会看到如何准备并启动直方图核函数。我们暂时假设已经在 GPU 上计算好了直方图。在计算完成后，需要将直方图复制回 CPU，因此我们分配了一个包含 256 个元素的数组，并且执行从设备到主机的复制。

unsigned int histo[256];
HANDLE_ERROR(cudaMemcpy(histo, dev_histo, 2

此时，我们完成了直方图的计算，因此可以停止计时器并显示经历的时间。

// 得到停止时间并显示计时结果
HANDLE_ERROR(cudaEventRecord(stop, 0));
HANDLE_ERROR(cudaEventSynchronize(stop));
float elapsedTime;
HANDLE_ERROR(cudaEventElapsedTime(&elapsedTime, start, stop));
std::cout << "Time to generate: " << elapsedTime << " ms\n";

此时，我们可以将直方图作为输入数据传递给算法的下一个步骤。然而，在这个示例中不需要将直方图用于其他任何操作，而只是验证在 GPU 上计算得到的直方图与在 CPU 上计算得到的直方图是否相同。首先，我们验证直方图的总和等于正确的值。这与 CPU 版本中的代码是相同的，如下所示：

long histoCount = 0;
for (int i = 0; i < 256; i++) {
     histoCount += histo[i];
}
std::cout << "Histogram Sum: " << histoCount << std::endl;

计算出 GPU 直方图，并在遍历每个数值时，递减直方图中相应元素的值。因此，如果完成计算时直方图每个元素的值都为 0，那么 CPU 计算的直方图与 GPU 计算的直方图相等。从某种意义上来说，我们是在计算这两个直方图之间的差异。

// 验证与基于CPU计算得到的结果是相同的
for (int i = 0; i < SIZE; i++)
     histo[buffer[i]]--;
for (int i = 0; i < 256; i++) {
    if (histo[i] != 0)
        std::cout << "Failure at " << i << "!\n";
}

程序结束前要施放已分配的 CUDA 事件，GPU 内存和主机内存。

	HANDLE_ERROR(cudaEventDestroy(start));
	HANDLE_ERROR(cudaEventDestroy(stop));
	cudaFree(dev_histo);
	cudaFree(dev_buffer);
	free(buffer);
	return 0;
}

由于直方图包含了 256 个元素，因此可以在每个线程块中包含 256 个线程，这种方式不仅方便而且高效。但是，在线程块的数量上还可以有更多选择。例如，在 100MB 数据中共有 104857600 个字节。我们可以启动一个线程块，并且让每个线程处理 409600 个数据元素。同样，我们还可以启动 409600 个线程块，并且让每个线程处理一个数据元素。

最优的解决方案是在这两种极端情况之间。通过一些性能实验，我们发现当线程块的数量为 GPU 中处理器数量的 2 倍时，将达到最优性能。例如，在 GeForce GTX280 中包含了 30 个处理器，因此当启动 60 个并行线程块时，直方图核函数将运行得最快。

如果要基于当前的硬件平台来动态调整线程块的数量，那么就要用到其中一个设备属性。我们通过以下代码片段来实现这个操作。

cudaDeviceProp prop;
HANDLE_ERROR(cudaGetDeviceProperties(&prop, 0));
int blocks = prop.multiProcessorCount;
std::cout << "MultiProcessorCount: " << blocks << std::endl;
histo_kernel << <blocks * 2, 256 >> > (dev_buffer, SIZE, dev_histo);

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

计算直方图的核函数需要的参数包括：

一个指向输入数组的指针
输入数组的长度
一个指向输出直方图的指针

核函数执行的第一个计算就是计算输入数据数组中的偏移。每个线程的起始偏移都是 0 到线程数量减 1 之间的某个值。然后，对偏移的增量为已启动线程的总数。


#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "../../common/book.h"
#include <stdio.h>
#include <iostream>
#define SIZE (100*1024*1024)

 // 使用全局内存原子操作的直方图核函数
__global__ void histo_kernel(unsigned char* buffer, long size, unsigned int* histo) {
    int i = threadIdx.x + blockIdx.x * blockDim.x;
    int stride = blockDim.x * gridDim.x;

    while (i < size) {
        atomicAdd(&(histo[buffer[i]]), 1);
        i += stride;
    }
}

函数调用 atomicAdd(addr, y) 将生成一个原子的操作序列，这个操作序列包括读取地址 addr 处的值，将 y 增加到这个值，以及将结果保存回地址 addr。

底层硬件将确保当前执行这些操作时，其他任何线程都不会读取或写入地址 addr 上的值，这样就能确保得到预计的结果。

在这里的示例中，这个地址就是直方图中相应元素的位置。如果当前字节为 buffer[i]，那么直方图中相应的元素就是 histo[buffer[i]]。原子操作需要这个元素的地址，因此第一个参数为 &(histo[buffer[i])。由于我们只是想把这个元素中的值递增 1，因此第二个参数就是 1。

完整代码


#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "../../common/book.h"
#include <stdio.h>
#include <iostream>
#define SIZE (100*1024*1024)

 //使用全局内存原子操作的直方图核函数
__global__ void histo_kernel(unsigned char* buffer, long size, unsigned int* histo) {
    int i = threadIdx.x + blockIdx.x * blockDim.x;
    int stride = blockDim.x * gridDim.x;

    while (i < size) {
        atomicAdd(&(histo[buffer[i]]), 1);
        i += stride;
    }
}

int main( void ) {
    unsigned char* buffer = (unsigned char*)big_random_block(SIZE);
    cudaEvent_t start, stop;
    HANDLE_ERROR(cudaEventCreate(&start));
    HANDLE_ERROR(cudaEventCreate(&stop));
    HANDLE_ERROR(cudaEventRecord(start, 0));

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

    cudaDeviceProp prop;
    HANDLE_ERROR(cudaGetDeviceProperties(&prop, 0));
    int blocks = prop.multiProcessorCount;
    std::cout << "MultiProcessorCount: " << blocks << std::endl;
    histo_kernel << <blocks * 2, 256 >> > (dev_buffer, SIZE, dev_histo);

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

    // 得到停止时间并显示计时结果
    HANDLE_ERROR(cudaEventRecord(stop, 0));
    HANDLE_ERROR(cudaEventSynchronize(stop));
    float elapsedTime;
    HANDLE_ERROR(cudaEventElapsedTime(&elapsedTime, start, stop));
    std::cout << "Time to generate: " << elapsedTime << " ms\n";

    long histoCount = 0;
    for (int i = 0; i < 256; i++) {
        histoCount += histo[i];
    }
    std::cout << "Histogram Sum: " << histoCount << std::endl;

    // 验证与基于CPU计算得到的结果是相同的
    for (int i = 0; i < SIZE; i++)
        histo[buffer[i]]--;
    for (int i = 0; i < 256; i++) {
        if (histo[i] != 0)
            std::cout << "Failure at " << i << "!\n";
    }

    HANDLE_ERROR(cudaEventDestroy(start));
    HANDLE_ERROR(cudaEventDestroy(stop));
    cudaFree(dev_histo);
    cudaFree(dev_buffer);
    free(buffer);
    return 0;
}

运行结果

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

由于在核函数中只包含了非少的计算工作，因此很可能全局内存上的原子操作导致性能的降低。

当数千个线程尝试访问少量的内存位置时，将发生大量的竞争。为了确保递增操作的原子性，对相同内存位置的操作都将被硬件串行化。这可能导致未完成操作的队列非常长，因此会抵消通过并行运行线程而获得的性能提升。

尽管这些原子操作是导致这种性能降低的原因，但解决这个问题的方法却出乎意料地需要使用更多而非更少的原子操作。

这里的主要问题并非在于使用了过多的原子操作，而是有数千个线程在少量的内存地址上发生竞争。要解决这个问题，我们将直方图计算分为两个阶段。

在第一个阶段，每个并行线程块将计算它所处理数据的直方图。由于每个线程块在执行这个操作时都是相互独立的，因此可以在共享内存中计算这些直方图，这将避免每次将写入操作从芯片发送到DRAM。但是，这种方式仍然需要原子操作，因为在线程块中多个线程之间仍然会处理相同值的数据元素。然而，现在只有 256 个线程在 256 个地址上发生竞争，这将极大地减少在使用全局内存时在数千个线程之间发生竞争的情况。

然后，在第一个阶段中分配一个共享内存缓冲区并进行初始化，用来保存每个线程块的临时直方图。由于随后的步骤将包括读取和修改这个缓冲区，因此需要调用 __syncthreads() 来确保每个线程的写入操作在线程继续前进之前完成。

__global__ void histo_kernel(unsigned char* buffer, long size, unsigned int* histo) {
    __shared__ unsigned int temp[256];
    temp[threadIdx.x] = 0;
    __syncthreads();

在将直方图初始化为 0 后，下一个步骤与最初 GPU 版本的直方图计算非常类似。这里唯一的差异在于，我们使用了共享内存缓冲区 temp[] 而不是全局内存缓冲区 histo[]，并且需要随后调用 __syncthreads() 来确保提交最后的写入操作。

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

    __syncthreads();

最后一步要求将每个线程块的临时直方图合并到全局缓冲区 histo[] 中。假设将输入数据分为两半，这样就有两个线程查看不同部分的数据，并计算得到两个独立的直方图。如果线程 A 在输入数据中发现字节 0xFC 出现了 20 次，线程 B 发现字节 0xFC 出现了 5 次，那么字节 0xFC 在输入数据中共出现了 25 次。同样，最终直方图的每个元素只是线程 A 直方图中相应元素和线程 B 直方图中相应元素的加和。

这个逻辑可以扩展到任意数量的线程，因此将每个线程块的直方图合并为单个最终的直方图就是，将线程块的直方图的每个元素都相加到最终直方图中相应位置的元素上。这个操作需要自动完成：

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

由于我们使用了 256 个线程，并且直方图中包含了 256 个元素，因此每个线程都将自动把它计算得到的元素只增加到最终直方图的元素上。如果线程数量不等于元素数量，那么这个阶段将更为复杂。

注意，我们并不保证线程块将按照何种顺序将各自的值相加到最终直方图中，但由于整数加法是可交换的，无论哪种顺序都会得到相同的结果。

完整代码


#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "../../common/book.h"
#include <stdio.h>
#include <iostream>
#define SIZE (100*1024*1024)

// 使用共享内存原子操作和全局内存原子操作的直方图核函数
__global__ void histo_kernel(unsigned char* buffer, long size, unsigned int* histo) {
    __shared__ unsigned int temp[256];
    temp[threadIdx.x] = 0;
    __syncthreads();

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

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

int main( void ) {
    unsigned char* buffer = (unsigned char*)big_random_block(SIZE);
    cudaEvent_t start, stop;
    HANDLE_ERROR(cudaEventCreate(&start));
    HANDLE_ERROR(cudaEventCreate(&stop));
    HANDLE_ERROR(cudaEventRecord(start, 0));

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

    cudaDeviceProp prop;
    HANDLE_ERROR(cudaGetDeviceProperties(&prop, 0));
    int blocks = prop.multiProcessorCount;
    std::cout << "MultiProcessorCount: " << blocks << std::endl;
    histo_kernel << <blocks * 2, 256 >> > (dev_buffer, SIZE, dev_histo);

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

    // 得到停止时间并显示计时结果
    HANDLE_ERROR(cudaEventRecord(stop, 0));
    HANDLE_ERROR(cudaEventSynchronize(stop));
    float elapsedTime;
    HANDLE_ERROR(cudaEventElapsedTime(&elapsedTime, start, stop));
    std::cout << "Time to generate: " << elapsedTime << " ms\n";

    long histoCount = 0;
    for (int i = 0; i < 256; i++) {
        histoCount += histo[i];
    }
    std::cout << "Histogram Sum: " << histoCount << std::endl;

    // 验证与基于CPU计算得到的结果是相同的
    for (int i = 0; i < SIZE; i++)
        histo[buffer[i]]--;
    for (int i = 0; i < 256; i++) {
        if (histo[i] != 0)
            std::cout << "Failure at " << i << "!\n";
    }

    HANDLE_ERROR(cudaEventDestroy(start));
    HANDLE_ERROR(cudaEventDestroy(stop));
    cudaFree(dev_histo);
    cudaFree(dev_buffer);
    free(buffer);
    return 0;
}