减少后的不同结果

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 

有很多事情要提。

  1. 我不确定您的错误检查是否有效。你没有显示实现这个的文件,当我 运行 你的代码 cuda-memcheck 时,我收到各种错误报告。好像都跟锁定功能有关

  2. 我不知道你为什么要使用锁定功能,我不推荐它。根据您使用它的方式,我认为它不会在您认为超出范围的情况下超出范围。我建议改用 atomicAdd ,这应该更快更简单。至少,在析构函数中注释掉 cudaFree() 语句。

  3. 您正在链接到旧的演示文稿。如果您查看 newer version of it,那么我想您会发现它现在建议使用 volatile。我不会为您重写您的整个代码,也不会再次总结该白皮书,但如果您只是为了演示目的将 volatile 添加到共享内存声明,它将解决由此产生的问题。

  4. 您的共享内存变量声明为 int,但您正在对 float 个数量求和。那不会按照你想要的方式工作。您可以这样声明它:

    extern __shared__ volatile float sdata[];
    
  5. 以上更改为我获取代码 "working"。剩下的一项是 CPU 和 GPU 结果之间的差异。我认为这是由于 CPU(串行缩减)与 GPU(并行缩减)的浮点运算顺序不同。由于差异出现在 float 数量的第 6 位有效数字,我认为这完全在浮点结果比较的合理范围内。如果您需要更高的准确性,您可以尝试从 float 切换到 double。此外,您可能希望阅读各种浮点论文,这将有助于理解此处,例如 this one and this one.