CUDA原子操作

考虑当大量的线程需要试图修改一段较小的内存区域的情形，这是(在日常的算法实现中）常发生的现象。当我们试图进行“读取–修改-写入”操作序列的时候，这种情形经常会带来很多麻烦。一个例子是代码d_out[i]++，这代码首先将d_out[i]的原值从存储器中读取出来，然后执行了+1操作，再将结果回写到存储器。然而，如果多个线程试图在同一个内存区域中进行这个操作，则可能会得到错误的结果。

#include <stdio.h>

#define NUM_THREADS 10000
#define SIZE  10
#define BLOCK_WIDTH 100

__global__ void gpu_increment_without_atomic(int *d_a)
{
	// Calculate thread id for current thread
	int tid = blockIdx.x * blockDim.x + threadIdx.x;

	// each thread increments elements wrapping at SIZE variable
	tid = tid % SIZE;
	d_a[tid] += 1;
}

int main()
{
	printf("%d total threads in %d blocks writing into %d array elements\n",
		NUM_THREADS, NUM_THREADS / BLOCK_WIDTH, SIZE);

	// declare and allocate host memory
	int h_a[SIZE];
	const int ARRAY_BYTES = SIZE * sizeof(int);

	// declare and allocate GPU memory
	int * d_a;
	cudaMalloc((void **)&d_a, ARRAY_BYTES);
	//Initialize GPU memory to zero
	cudaMemset((void *)d_a, 0, ARRAY_BYTES);

	gpu_increment_without_atomic << <NUM_THREADS / BLOCK_WIDTH, BLOCK_WIDTH >> >(d_a);

	// copy back the array to host memory
	cudaMemcpy(h_a, d_a, ARRAY_BYTES, cudaMemcpyDeviceToHost);

	printf("Number of times a particular Array index has been incremented without atomic add is: \n");
	for (int i = 0; i < SIZE; i++)
	{
		printf("index: %d --> %d times\n ", i, h_a[i]);
	}

	cudaFree(d_a);
	return 0;
}

内核函数简单地通过d_a[tid]+=1这行代码来增加存储器中元素的值。关键的问题在于，(这行代码）对应的具体内存区域被增加了多少次?线程总数为10 000，数组里只有10个(元素)。通过求余(求模，%)运算，来将这10 000个线程ID对应的索引到这10个元素上去。所以，每个相同的内存中的元素位置将有1 000个线程来进行（+1)的运算。理想状态下，数组中每个位置的元素将被增加1000(次个1)。但是我们从输出结果中看并非如此。
为了解决这个问题，CUDA提供了atomicAdd这种原子操作函数。该函数会从逻辑上保证，每个调用它的线程对相同的内存区域上的“读取旧值-累加-回写新值”操作是不可被其他线程扰乱的原子性的整体完成的。

#include <stdio.h>
#include <cuda.h>
#include <cuda_runtime.h>

#define NUM_THREADS 10000
#define SIZE  10
#define BLOCK_WIDTH 100

__global__ void gpu_increment_atomic(int *d_a)
{
	// Calculate thread id for current thread
	int tid = blockIdx.x * blockDim.x + threadIdx.x;

	// each thread increments elements wrapping at SIZE variable
	tid = tid % SIZE;
	atomicAdd(&d_a[tid], 1);
}

int main()
{
	printf("%d total threads in %d blocks writing into %d array elements\n",
		NUM_THREADS, NUM_THREADS / BLOCK_WIDTH, SIZE);

	// declare and allocate host memory
	int h_a[SIZE];
	const int ARRAY_BYTES = SIZE * sizeof(int);

	// declare and allocate GPU memory
	int * d_a;
	cudaMalloc((void **)&d_a, ARRAY_BYTES);
	//Initialize GPU memory to zero
	cudaMemset((void *)d_a, 0, ARRAY_BYTES);
	
	gpu_increment_atomic << <NUM_THREADS / BLOCK_WIDTH, BLOCK_WIDTH >> >(d_a);
	
	// copy back the array to host memory
	cudaMemcpy(h_a, d_a, ARRAY_BYTES, cudaMemcpyDeviceToHost);
	
	printf("Number of times a particular Array index has been incremented is: \n");
	for (int i = 0; i < SIZE; i++) 
	{ 
		printf("index: %d --> %d times\n ", i, h_a[i]); 
	}
	
	cudaFree(d_a);
	return 0;
}

我们用atomicAdd原子操作函数替换了之前的直接+=操作，该函数具有2个参数:第一个参数是我们要进行原子加法操作的内存区域;第二个参数是该原子加法操作具体要加上的值。在这个代码中，1 000个线程对同一内存区域进行原子+1操作，这1 000次相同区域上的操作，每次都将从逻辑上安全地完整执行。这可能会增加执行时间上的代价。
在main函数中，具有10个元素的数组被初始化成0值，然后传递给了内核，但现在，内核中的代码将执行原子累加操作。所以，这个程序输出的结果将是对的，数组中的每个元素将被累加1000。
如果你测量一下这个程序的运行时间，相比之前的那个简单地在全局内存上直接进行加法操作的程序它用的时间更长。这是因为使用原子操作后程序具有更大的执行代价。可以通过使用共享内存来加速这些原子累加操作。如果线程规模不变，但原子操作的元素数量扩大，则这些同样次数的原子操作会更快地完成。这是因为更广泛的分布范围上的原子操作有利于利用多个能执行原子操作的单元，以及每个原子操作单元上面的竞争性的原子事务也相应减少了。