CUDA By Example(七)——原子性

阿里云国内75折回扣微信号：monov8

阿里云国际，腾讯云国际，低至75折。AWS 93折免费开户实名账号代冲值优惠多多微信号：monov8 飞机：@monov6

文章目录

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

原子操作简介

在编写传统的单线程应用程序时程序员通常不需要使用原子操作。即使需要使用原子操作也无需担心。下面将详细解释原子操作是什么以及为什么在多线程程序中需要使用它们。为了解释原子操作首先来分析 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;
}

运行结果

可以看到相比于仅仅使用全局内存原子操作的版本要快了 1 倍左右。

CUDA By Example(七)——原子性

文章目录

原子操作简介

计算直方图

在 CPU 上计算直方图

在GPU上计算直方图

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

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

更多相关内容

“CUDA By Example(七)——原子性” 的相关文章

阿里云国际版