CUDA计算直方图，核函数启动二维线程块、线程格，但是只能计算size大小为16的倍数的data

[dreamcraneleaf]

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

#define M 10240
#define N 10240
#define SIZE (M*N)

// Simple utility function to check for CUDA runtime errors
void checkCUDAError(const char *msg)
{
cudaError_t err = cudaGetLastError();
if( cudaSuccess != err)
{
fprintf(stderr, "Cuda error: %s: %s.\n", msg, cudaGetErrorString( err) );
getchar();
exit(-1);
}
}

__global__ void myhistKernel(unsigned char * buffer,unsigned int * histo)
{
__shared__ unsigned int temp[256];

int x = blockIdx.x * blockDim.x + threadIdx.x;
int y = blockIdx.y * blockDim.y + threadIdx.y;
int offset = x + y * blockDim.x * gridDim.x;
int linear = threadIdx.x + threadIdx.y * blockDim.x;
if(linear<256)
temp[linear] = 0;
__syncthreads();

if(x<M && y<N)
atomicAdd( &temp[buffer[offset]], 1 );

__syncthreads();
if(linear<256)
atomicAdd(&histo[linear], temp[linear]);
}

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 ) ) );

dim3 dimBlock(16,16);
dim3 dimGrid((10240+dimBlock.x-1)/(dimBlock.x),(10240+dimBlock.y-1)/(dimBlock.y));
myhistKernel<<<dimGrid,dimBlock>>>(dev_buffer,dev_histo);
cudaThreadSynchronize();
checkCUDAError("kernel execution");

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 );
getchar();
return 0;
}
如果data不是16的倍数，buffer[offset]将出现内存异常，哪位朋友知道采取什么样的方法可以解决这个问题，谢谢

[Spidey212]

dim3 dimGrid((10240+dimBlock.x-1)/(dimBlock.x),(10240+dimBlock.y-1)/(dimBlock.y)); 这一句改成

dim3 dimGrid((M+dimBlock.x-1)/(dimBlock.x),(N+dimBlock.y-1)/(dimBlock.y));

太粗心了

[dreamcraneleaf]

这样写没解决问题啊，大神！如果M=N=1000，不是16的倍数，这段代码是要出问题的，主要原因是buffer[offset]获取的值为空

[Spidey212]

实际上，当M不为16的整数倍时会出错，而N则没有这个限制。因为

 int offset = x + y * blockDim.x * gridDim.x; 

去buffer里取数据的线程网格尺寸是 （blockDim.x * gridDim.x）* N； 后面的if语句

	
if(x<M && y<N)
		atomicAdd( &temp[buffer[offset]], 1 );

限制的二维线程块是M*N，这两个显然不匹配。 改成

int offset = x + y * M;

[dreamcraneleaf]

非常感谢您的帮助，问题终于解决了，但是执行效率比较低，想咨询一下我怎么样做才可以提高计算效率？个人在尝试分配显存时候采用cudaMallocPitch()和cudaMemcpy2D(),不过总是提示内存冲突，还望得到您的指点。非常感谢，祝好

[dreamcraneleaf]

非常感谢您这几天的帮助，打扰您了。接下来我会继续努力去优化基于GPU的直方图计算，你给的建议我也会去研究，以后会多向您学习的。

[Spidey212]

引用 4 楼 dreamcraneleaf 的回复:

非常感谢您的帮助，问题终于解决了，但是执行效率比较低，想咨询一下我怎么样做才可以提高计算效率？个人在尝试分配显存时候采用cudaMallocPitch()和cudaMemcpy2D(),不过总是提示内存冲突，还望得到您的指点。非常感谢，祝好

原子操作效率很低，当多个线程发生数据冲突时会用原子操作来让线程互斥，原子函数的操作是串行的，尽量少用。 代码里第一个还好，毕竟随机数分布比较平均，同时写一个地址的线程不会太多。 而第二个atomicAdd(&histo[linear], temp[linear]);一共有多少个block，就会有多少线程写同一地址，当数据规模较小时，这样做没问题，这里可以估算一下，每个block有256个线程，这里面肯定有256个原子操作并行进行的，发生数据冲突的线程数小于等于256时，效率不会受到太大影响，所以block的维度不能高于16 * 16，也就得出数据规模的阈值是16 * 256 = 4096。当M = N <= 4096时，速度应该是可以接受的，再往上的话，数据规模越大，效率会越来越低，就需要考虑优化了。 你可以实际测试一下。 优化方案就可以考虑取消原方案的原子操作，用归约来替代；甚至，可以应用MapReduce的思想来做，自己找找相关资料。 PS：优化方案可能会带来任务调度，内存分配以及kernel执行的额外时间开销，所以实际应用中应该根据数据规模来选择合适的方案