简介
第9章引入了3个新的概念，在之前的章节和示例中都尽量避免这方面的内容。
- 首先是Compute Capability. NVIDIA不断发布新的GPU，计算能力增大的同时，还伴随有新功能的添加。比如本章中要介绍的原子操作。Compute Capability 1.0不支持原子操作，1.1支持Global Memory原子操作，1.2支持Shared Momory原子操作。所以要使用某项功能时，一定要对运行的硬件环境心中有数。GPU型号的Compute Capability可以在NVIDIA的官网上查到。
- 其次是原子操作。我们对于CPU上的原子操作不会陌生。当多个线程对某个资源同时操作时，要确保一个线程完成操作之后，才让另外的线程来操作。原子操作可以保证做到这一点。前面的示例，都尽量让每个线程操作对应的内存位置，避免冲突。
- 直方图。直方图我们可以简单地理解为，将事物进行分类，然后统计各个分类中的数量。我们要确定好分成多少类，分类的标准是什么。本章的示例，是将生成的大量数据进行分类，数据长度为1个字节。然后将数据分成256类，每个数据的数值大小作为分类的标准。
直方图CPU示例
数据获取
直方图示例其实挺简单的。唯一麻烦一点的数据从哪里来。Cuda By Example 的随书代码实现了一个生成大块数据的函数big_random_block进一步简化我们的工作。
#include "../common/book.h"
#define SIZE (100*1024*1024)
int main( void ) {
unsigned char *buffer = (unsigned char*)big_random_block( SIZE );生成直方图
然后定义一个长度为256的无符号整型数组，对应直方图的256个区间，统计0~255在buffer出现次数。
unsigned int histo[256];
for (int i=0; i<256; i++)
histo[i] = 0;
for (int i=0; i<SIZE; i++)
histo[buffer[i]]++;histo数组的各个元素被第一个for循环初始化为0，第二个for循环遍历buffer中的元素，buffer[i]的计数放在histo[buffer[i]]中。
验证结果
生成直方图的过程到此就结束了，再做一个简单的验证，将各个区间的计数求和，所得结果应该跟buffer的SIZE相同。
long histoCount = 0;
for (int i = 0; i<256; i++)
histoCount += histo[i];
printf("Histogram Sum: %ld\n", histoCount);测量性能
GPU的性能我们使用CUDA提供的CudaEvent来测量， Host端的性能测量要简单得多。标准C语言库提供了一个函数clock()，通过调用此函数，我们可以取得系统自开机以来，经历的时钟周期数。在不同时刻调用这个函数，所得返回值的差即为两个时刻相隔的时钟周期数量。时钟周期数量除以每秒时钟数，我们就得到了以秒为单位的时间间隔。生成直方图消耗的时间不会太多，所以取毫秒为单位，这样还需将所得结果乘以1000。
clock()返回的数据类型clock_t实际上由整型或者浮点数的typedef而来，可以直接使用减号来计算差值。
#include <timer.h>
clock_t start, stop;
start = clock();
...
stop = clock();
float elapsedTime = (float)(stop-start)/(float)CLOCKS_PER_SEC*1000.0f;
printf( "Time to generate: %3.1f ms\n", elapsedTime );完整代码
#include "../common/book.h"
#define SIZE (100*1024*1024)
int main( void ) {
unsigned char *buffer =
(unsigned char*)big_random_block( SIZE );
// capture the start time
clock_t start, stop;
start = clock();
unsigned int histo[256];
for (int i=0; i<256; i++)
histo[i] = 0;
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 );
long histoCount = 0;
for (int i=0; i<256; i++) {
histoCount += histo[i];
}
printf( "Histogram Sum: %ld\n", histoCount );
free( buffer );
return 0;
}
直方图GPU示例（使用全局变量）
内核函数
使用GPU生成直方图，让不同的GPU线程，检查buffer的不同部分。这个前面的例子已经讲过很多次了。这个例子不同的地方是，把结果写回histo时，线程间可能同时访问某个histo区间，从而产生冲突。为了避免出现统计错误，线程必须使用原子操作来改变histo各个元素的值。CUDA C用于加法的原子操作为：atomicAdd(addr, y)。addr为要进行原子操作的变量的地址，y为要加入到变量的数值。例子中addr为histo里buffer[i] 的计数的地址，&histo[buffer[i]]。y为1，因为对于数值buffer[i]，计数histo[buffer[i]]加1.
__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;
}
}数据准备
正如之前例子所作的那样，需要分配GPU内存保存buffer数据，并将数据从Host内存拷贝到GPU内存。同时也需要给histogram分配内存，供内核线程存放统计数据。这里不要忘了将histo初始化为0. CUDA C的cudaMemset()起到跟标准C的memset类似的作用。例子使用cudaMemset初始化了histo。
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( hist_histo, 0, 256 * sizeof(int) ) );返回histo
初始化输入和输出内存后，就可以进行直方图的计算了。就像用malloc分配了内存后，马上把释放的代码也写好是个好习惯，这里在写出计算代码之前，我们把从GPU读回结果先介绍了。GPU计算的直方图结果将传递给Host的临时数组histo，用于后面对计算结果的验证。
unsigned int histo[256];
HANDLE_ERROR( cudaMemcpy( histo, dev_histo,
256 * sizeof( int ),
cudaMemcpyDeviceToHost ) );验证结果
CPU版本里， 验证过程是计算各个区间的和，然后把和跟数据总数量做比较，如果相同则认为是准确的。GPU版本也可以使用同样的方法来验证，不过这里有更好的方法。以GPU计算所得的直方图为基础，在Host端重新遍历输入buffer，每次从对应的区间中减去1. 遍历完之后，histo的每个元素应该全部归为零值。
for (int i=0; i<SIZE; i++)
histo[buffer[i]]--;
for (int i=0; i<256; i++) {
if (histo[i]) != 0)
printf("Failare at %d!\n", i);
} 测量性能
跟之前一样，例子使用start和stop两个cuda Event来测量直方图计算所花费的时间。因为跟前面的例子完全一样，这里不再赘述了。具体的代码可以在后面的完整代码中找到。
确定blocks
将调用kernel函数放到最后面将是有原因的，因为选择不同的blocks对性能是由影响的。除了我们首先想到的每个block使用256个thread，从而计算得出block数量，我们还有其他的选择。实验证明，将block数量定义为GPU包含的处理器的2倍，性能是最好的。
我们还是使用每个block建立256个线程，block由GPU的prop来决定。cudaGetDeviceProperties函数返回GPU的属性，用cduaDeviceProp来描述。其中multiProcessorCount指定了GPU里处理器的数量。这样就得到了下面的调用代码：
cudaDeviceProp prop;
HANDLE_ERROR( cudaGetDeviceProperties( &prop, 0 ) );
int blocks = prop.multiProcessorCount;
histo_kernel<<<blocks*2, 256>>>(dev_buffer, SIZE, dev_histo ));完整代码
#include "../common/book.h"
#define SIZE (100*1024*1024)
__global__ void histo_kernel( unsigned char *buffer,
long size,
unsigned int *histo ) {
// calculate the starting index and the offset to the next
// block that each thread will be processing
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 );
// capture the start time
// starting the timer here so that we include the cost of
// all of the operations on the GPU.
cudaEvent_t start, stop;
HANDLE_ERROR( cudaEventCreate( &start ) );
HANDLE_ERROR( cudaEventCreate( &stop ) );
HANDLE_ERROR( cudaEventRecord( start, 0 ) );
// allocate memory on the GPU for the file's data
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 ) ) );
// kernel launch - 2x the number of mps gave best timing
cudaDeviceProp prop;
HANDLE_ERROR( cudaGetDeviceProperties( &prop, 0 ) );
int blocks = prop.multiProcessorCount;
histo_kernel<<<blocks*2,256>>>( dev_buffer, SIZE, dev_histo );
unsigned int histo[256];
HANDLE_ERROR( cudaMemcpy( histo, dev_histo,
256 * sizeof( int ),
cudaMemcpyDeviceToHost ) );
// get stop time, and display the timing results
HANDLE_ERROR( cudaEventRecord( stop, 0 ) );
HANDLE_ERROR( cudaEventSynchronize( stop ) );
float elapsedTime;
HANDLE_ERROR( cudaEventElapsedTime( &elapsedTime,
start, stop ) );
printf( "Time to generate: %3.1f ms\n", elapsedTime );
long histoCount = 0;
for (int i=0; i<256; i++) {
histoCount += histo[i];
}
printf( "Histogram Sum: %ld\n", histoCount );
// verify that we have the same counts via 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( stop ) );
cudaFree( dev_histo );
cudaFree( dev_buffer );
free( buffer );
return 0;
}性能测量得出的时间相比CPU版本可能不尽人意. 因为dev_histo是全局变量,计算的时候多个线程会访问dev_histo的同一个元素. 比如对于&dev_histo[1], 所有线程,无论属于哪个block, 如果它从dev_buffer中读到数据为1, 那么都要来排队访问&dev_histo[1].内存访问极大影响了性能.通过前面的学习,我们知道共享内存可以很好解决这一问题,因为共享内存不仅访问速度比全局内存快得多,而且只有同一个block的线程来访问, 比起所有的线程,竞争小多了. 只是最后需要把所有block的计数加在一起,不过这不算什么问题.
直方图GPU示例（使用共享内存）
内核函数
使用共享内存,只需要修改内核函数即可, 其他的地方不需要动. 不过需要注意, 每个block的线程数需要设定为256个,因为在下面的内核函数中, 共享内存的初始化和将计算所得结果加到dev_histo中去, 都是一个thread负责一个区间. 如果小于256则会漏掉一些区间, 如果大于256则会导致访问越界.
__global__ void histo_kernel( unsigned char* dev_buffer,
long size,
unsigned int *dev_histo)
{
__shared__ unsigned int temp[256];
temp[threadIdx.x] = 0;
__syncthreads();
int i = threadIdx.x + blockIdx.x * blockDim.x;
int stride = blockDim.x * gridDim.x;
while (i < size ) {
atomicAdd( &temp[dev_buffer[i]], 1);
i += stride;
}
__syncthreads();
atomicAdd(&dev_histo[threadIdx.x], temp[threadIdx.x]);
}
这个函数中间的部分,跟全局版本差不多,只是计数结果不在直接放到dev_histo中,而是放入定义在函数内部的共享内存中.
开头部分申请了共享内存 temp. 对共享内存的初始化, 每个线程负责一个位置, 由线程在block中索引作为数组的索引. 而之后的__syncthreads()则是为了保证每个线程都运行了temp[threadIdx.x]=0, 将整个数组的每个元素清零后, 才开始后面运算.
第二个__syncthreads()确保了block中所有线程都完成了计算, temp得到最终值后,才加到dev_histo中去.同样每个线程负责一个区间的叠加.
完整代码
虽然改动不多,但是还是贴上完整的代码,毕竟说得再多,都不如代码直接.
#include "../common/book.h"
#define SIZE (100*1024*1024)
__global__ void histo_kernel( unsigned char *buffer,
long size,
unsigned int *histo ) {
// clear out the accumulation buffer called temp
// since we are launched with 256 threads, it is easy
// to clear that memory with one write per thread
__shared__ unsigned int temp[256];
temp[threadIdx.x] = 0;
__syncthreads();
// calculate the starting index and the offset to the next
// block that each thread will be processing
int i = threadIdx.x + blockIdx.x * blockDim.x;
int stride = blockDim.x * gridDim.x;
while (i < size) {
atomicAdd( &temp[buffer[i]], 1 );
i += stride;
}
// sync the data from the above writes to shared memory
// then add the shared memory values to the values from
// the other thread blocks using global memory
// atomic adds
// same as before, since we have 256 threads, updating the
// global histogram is just one write per thread!
__syncthreads();
atomicAdd( &(histo[threadIdx.x]), temp[threadIdx.x] );
}
int main( void ) {
unsigned char *buffer =
(unsigned char*)big_random_block( SIZE );
// capture the start time
// starting the timer here so that we include the cost of
// all of the operations on the GPU. if the data were
// already on the GPU and we just timed the kernel
// the timing would drop from 74 ms to 15 ms. Very fast.
cudaEvent_t start, stop;
HANDLE_ERROR( cudaEventCreate( &start ) );
HANDLE_ERROR( cudaEventCreate( &stop ) );
HANDLE_ERROR( cudaEventRecord( start, 0 ) );
// allocate memory on the GPU for the file's data
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 ) ) );
// kernel launch - 2x the number of mps gave best timing
cudaDeviceProp prop;
HANDLE_ERROR( cudaGetDeviceProperties( &prop, 0 ) );
int blocks = prop.multiProcessorCount;
histo_kernel<<<blocks*2,256>>>( dev_buffer,
SIZE, dev_histo );
unsigned int histo[256];
HANDLE_ERROR( cudaMemcpy( histo, dev_histo,
256 * sizeof( int ),
cudaMemcpyDeviceToHost ) );
// get stop time, and display the timing results
HANDLE_ERROR( cudaEventRecord( stop, 0 ) );
HANDLE_ERROR( cudaEventSynchronize( stop ) );
float elapsedTime;
HANDLE_ERROR( cudaEventElapsedTime( &elapsedTime,
start, stop ) );
printf( "Time to generate: %3.1f ms\n", elapsedTime );
long histoCount = 0;
for (int i=0; i<256; i++) {
histoCount += histo[i];
}
printf( "Histogram Sum: %ld\n", histoCount );
// verify that we have the same counts via 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!\n", i );
}
HANDLE_ERROR( cudaEventDestroy( start ) );
HANDLE_ERROR( cudaEventDestroy( stop ) );
cudaFree( dev_histo );
cudaFree( dev_buffer );
free( buffer );
return 0;
}本站资源均来自互联网，仅供研究学习，禁止违法使用和商用，产生法律纠纷本站概不负责！如果侵犯了您的权益请与我们联系！
转载请注明出处： 免费源码网-免费的源码资源网站 » Cuda By Example - 13 (原子操作）
发表评论 取消回复