减少后的不同结果
Different results after reduction
我有两个缩减算法,都来自 docs.nvidia,所以它们应该是正确的,但第一个(非常有效)给了我一个错误的结果。第二个结果更好,但我期望更高的准确性。是算法有问题还是我做的不好?
#include <stdio.h>
#include <cuda.h>
#include <stdlib.h>
#include <math.h>
#include "cuda_error.h"
//Lock definition
#ifndef __LOCK_H__
#define __LOCK_H__
struct Lock {
int *mutex;
Lock( void ) {
CudaSafeCall( cudaMalloc( (void**)&mutex,
sizeof(int) ) );
CudaSafeCall( cudaMemset( mutex, 0, sizeof(int) ) );
}
~Lock( void ) {
cudaFree( mutex );
}
__device__ void lock( void ) {
while( atomicCAS( mutex, 0, 1 ) != 0 );
}
__device__ void unlock( void ) {
atomicExch( mutex, 0 );
}
};
#endif
//-------------------------
const int N = 507904;
const int blockSize = 256;
int blocks = N/blockSize;
template <unsigned int blockSize>
__global__ void reduce(Lock lock, float *g_idata, float *g_odata, int n)
{
extern __shared__ int sdata[];
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x*(blockSize*2) + tid;
unsigned int gridSize = blockSize*2*gridDim.x;
sdata[tid] = 0;
while (i < n)
{
sdata[tid] += g_idata[i] + g_idata[i+blockSize];
i += gridSize;
}
__syncthreads();
if (blockSize >= 512) { if (tid < 256) { sdata[tid] += sdata[tid + 256]; } __syncthreads(); }
if (blockSize >= 256) { if (tid < 128) { sdata[tid] += sdata[tid + 128]; } __syncthreads(); }
if (blockSize >= 128) { if (tid < 64) { sdata[tid] += sdata[tid + 64]; } __syncthreads(); }
if (tid < 32)
{
if (blockSize >= 64) sdata[tid] += sdata[tid + 32];
if (blockSize >= 32) sdata[tid] += sdata[tid + 16];
if (blockSize >= 16) sdata[tid] += sdata[tid + 8];
if (blockSize >= 8) sdata[tid] += sdata[tid + 4];
if (blockSize >= 4) sdata[tid] += sdata[tid + 2];
if (blockSize >= 2) sdata[tid] += sdata[tid + 1];
}
if (tid == 0)
{
lock.lock();
*g_odata += sdata[0];
lock.unlock();
}
}
__global__ void reduction_sum(Lock lock,float *in, float *out, int N)
{
extern __shared__ float sf_data[];
int tid = threadIdx.x;
int i = blockIdx.x * blockDim.x + threadIdx.x;
sf_data[tid] = (i<N) ? in[i] : 0;
__syncthreads();
for (int s = blockDim.x/2; s>0; s>>=1)
{
if (tid < s)
{
sf_data[tid] += sf_data[tid + s];
}
__syncthreads();
}
if (tid == 0)
{
lock.lock();
*out += sf_data[0];
lock.unlock();
}
}
//initializing function
double los()
{
return (double)rand()/(double)RAND_MAX;
}
//-------------------------------------------
int main (void)
{
//memory allocation
float *a;
float *dev_a, *dev_blocksum1, *dev_blocksum2;
float s1=0, s2=0, spr=0;
a = (float*)malloc( N*sizeof(float) );
CudaSafeCall( cudaMalloc( (void**)&dev_a, N*sizeof(float) ) );
CudaSafeCall( cudaMemset( dev_a, 0, N*sizeof(float) ) );
CudaSafeCall( cudaMalloc( (void**)&dev_blocksum1, sizeof(float) ) );
CudaSafeCall( cudaMalloc( (void**)&dev_blocksum2, sizeof(float) ) );
CudaSafeCall( cudaMemset( dev_blocksum1, 0, sizeof(float) ) );
CudaSafeCall( cudaMemset( dev_blocksum2, 0, sizeof(float) ) );
//--------------------
//drawing, array fill
srand(2403);
int i;
for (i=0; i<N; i++)
{
a[i]=los();
spr+=a[i];
}
printf("CPU sum: %f \n", spr);
//copy HtoD
CudaSafeCall( cudaMemcpy( dev_a, a, N*sizeof(float), cudaMemcpyHostToDevice ) );
//---------------------
Lock lock;
//reduce
reduce<blockSize><<<blocks, blockSize, blockSize*sizeof(float)>>>(lock, dev_a, dev_blocksum1, N);
CudaSafeCall( cudaMemcpy ( &s1, dev_blocksum1, sizeof(float), cudaMemcpyDeviceToHost ) );
printf("GPU sum1: %f \n", s1);
//-----------------------
//reduction_sum
reduction_sum<<<blocks, blockSize, blockSize*sizeof(float)>>>(lock, dev_a, dev_blocksum2, N);
CudaSafeCall( cudaMemcpy ( &s2, dev_blocksum2, sizeof(float), cudaMemcpyDeviceToHost ) );
printf("GPU sum2: %f \n", s2);
//---------------------
return 0;
}
$ CPU sum: 253833.515625
$ GPU sum1: 14021.000000
$ GPU sum2: 253835.906250
有很多事情要提。
我不确定您的错误检查是否有效。你没有显示实现这个的文件,当我 运行 你的代码 cuda-memcheck
时,我收到各种错误报告。好像都跟锁定功能有关
我不知道你为什么要使用锁定功能,我不推荐它。根据您使用它的方式,我认为它不会在您认为超出范围的情况下超出范围。我建议改用 atomicAdd
,这应该更快更简单。至少,在析构函数中注释掉 cudaFree()
语句。
您正在链接到旧的演示文稿。如果您查看 newer version of it,那么我想您会发现它现在建议使用 volatile
。我不会为您重写您的整个代码,也不会再次总结该白皮书,但如果您只是为了演示目的将 volatile
添加到共享内存声明,它将解决由此产生的问题。
您的共享内存变量声明为 int
,但您正在对 float
个数量求和。那不会按照你想要的方式工作。您可以这样声明它:
extern __shared__ volatile float sdata[];
以上更改为我获取代码 "working"。剩下的一项是 CPU 和 GPU 结果之间的差异。我认为这是由于 CPU(串行缩减)与 GPU(并行缩减)的浮点运算顺序不同。由于差异出现在 float
数量的第 6 位有效数字,我认为这完全在浮点结果比较的合理范围内。如果您需要更高的准确性,您可以尝试从 float
切换到 double
。此外,您可能希望阅读各种浮点论文,这将有助于理解此处,例如 this one and this one.
我有两个缩减算法,都来自 docs.nvidia,所以它们应该是正确的,但第一个(非常有效)给了我一个错误的结果。第二个结果更好,但我期望更高的准确性。是算法有问题还是我做的不好?
#include <stdio.h>
#include <cuda.h>
#include <stdlib.h>
#include <math.h>
#include "cuda_error.h"
//Lock definition
#ifndef __LOCK_H__
#define __LOCK_H__
struct Lock {
int *mutex;
Lock( void ) {
CudaSafeCall( cudaMalloc( (void**)&mutex,
sizeof(int) ) );
CudaSafeCall( cudaMemset( mutex, 0, sizeof(int) ) );
}
~Lock( void ) {
cudaFree( mutex );
}
__device__ void lock( void ) {
while( atomicCAS( mutex, 0, 1 ) != 0 );
}
__device__ void unlock( void ) {
atomicExch( mutex, 0 );
}
};
#endif
//-------------------------
const int N = 507904;
const int blockSize = 256;
int blocks = N/blockSize;
template <unsigned int blockSize>
__global__ void reduce(Lock lock, float *g_idata, float *g_odata, int n)
{
extern __shared__ int sdata[];
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x*(blockSize*2) + tid;
unsigned int gridSize = blockSize*2*gridDim.x;
sdata[tid] = 0;
while (i < n)
{
sdata[tid] += g_idata[i] + g_idata[i+blockSize];
i += gridSize;
}
__syncthreads();
if (blockSize >= 512) { if (tid < 256) { sdata[tid] += sdata[tid + 256]; } __syncthreads(); }
if (blockSize >= 256) { if (tid < 128) { sdata[tid] += sdata[tid + 128]; } __syncthreads(); }
if (blockSize >= 128) { if (tid < 64) { sdata[tid] += sdata[tid + 64]; } __syncthreads(); }
if (tid < 32)
{
if (blockSize >= 64) sdata[tid] += sdata[tid + 32];
if (blockSize >= 32) sdata[tid] += sdata[tid + 16];
if (blockSize >= 16) sdata[tid] += sdata[tid + 8];
if (blockSize >= 8) sdata[tid] += sdata[tid + 4];
if (blockSize >= 4) sdata[tid] += sdata[tid + 2];
if (blockSize >= 2) sdata[tid] += sdata[tid + 1];
}
if (tid == 0)
{
lock.lock();
*g_odata += sdata[0];
lock.unlock();
}
}
__global__ void reduction_sum(Lock lock,float *in, float *out, int N)
{
extern __shared__ float sf_data[];
int tid = threadIdx.x;
int i = blockIdx.x * blockDim.x + threadIdx.x;
sf_data[tid] = (i<N) ? in[i] : 0;
__syncthreads();
for (int s = blockDim.x/2; s>0; s>>=1)
{
if (tid < s)
{
sf_data[tid] += sf_data[tid + s];
}
__syncthreads();
}
if (tid == 0)
{
lock.lock();
*out += sf_data[0];
lock.unlock();
}
}
//initializing function
double los()
{
return (double)rand()/(double)RAND_MAX;
}
//-------------------------------------------
int main (void)
{
//memory allocation
float *a;
float *dev_a, *dev_blocksum1, *dev_blocksum2;
float s1=0, s2=0, spr=0;
a = (float*)malloc( N*sizeof(float) );
CudaSafeCall( cudaMalloc( (void**)&dev_a, N*sizeof(float) ) );
CudaSafeCall( cudaMemset( dev_a, 0, N*sizeof(float) ) );
CudaSafeCall( cudaMalloc( (void**)&dev_blocksum1, sizeof(float) ) );
CudaSafeCall( cudaMalloc( (void**)&dev_blocksum2, sizeof(float) ) );
CudaSafeCall( cudaMemset( dev_blocksum1, 0, sizeof(float) ) );
CudaSafeCall( cudaMemset( dev_blocksum2, 0, sizeof(float) ) );
//--------------------
//drawing, array fill
srand(2403);
int i;
for (i=0; i<N; i++)
{
a[i]=los();
spr+=a[i];
}
printf("CPU sum: %f \n", spr);
//copy HtoD
CudaSafeCall( cudaMemcpy( dev_a, a, N*sizeof(float), cudaMemcpyHostToDevice ) );
//---------------------
Lock lock;
//reduce
reduce<blockSize><<<blocks, blockSize, blockSize*sizeof(float)>>>(lock, dev_a, dev_blocksum1, N);
CudaSafeCall( cudaMemcpy ( &s1, dev_blocksum1, sizeof(float), cudaMemcpyDeviceToHost ) );
printf("GPU sum1: %f \n", s1);
//-----------------------
//reduction_sum
reduction_sum<<<blocks, blockSize, blockSize*sizeof(float)>>>(lock, dev_a, dev_blocksum2, N);
CudaSafeCall( cudaMemcpy ( &s2, dev_blocksum2, sizeof(float), cudaMemcpyDeviceToHost ) );
printf("GPU sum2: %f \n", s2);
//---------------------
return 0;
}
$ CPU sum: 253833.515625
$ GPU sum1: 14021.000000
$ GPU sum2: 253835.906250
有很多事情要提。
我不确定您的错误检查是否有效。你没有显示实现这个的文件,当我 运行 你的代码
cuda-memcheck
时,我收到各种错误报告。好像都跟锁定功能有关我不知道你为什么要使用锁定功能,我不推荐它。根据您使用它的方式,我认为它不会在您认为超出范围的情况下超出范围。我建议改用
atomicAdd
,这应该更快更简单。至少,在析构函数中注释掉cudaFree()
语句。您正在链接到旧的演示文稿。如果您查看 newer version of it,那么我想您会发现它现在建议使用
volatile
。我不会为您重写您的整个代码,也不会再次总结该白皮书,但如果您只是为了演示目的将volatile
添加到共享内存声明,它将解决由此产生的问题。您的共享内存变量声明为
int
,但您正在对float
个数量求和。那不会按照你想要的方式工作。您可以这样声明它:extern __shared__ volatile float sdata[];
以上更改为我获取代码 "working"。剩下的一项是 CPU 和 GPU 结果之间的差异。我认为这是由于 CPU(串行缩减)与 GPU(并行缩减)的浮点运算顺序不同。由于差异出现在
float
数量的第 6 位有效数字,我认为这完全在浮点结果比较的合理范围内。如果您需要更高的准确性,您可以尝试从float
切换到double
。此外,您可能希望阅读各种浮点论文,这将有助于理解此处,例如 this one and this one.