带有 CUDA 内核的点积用于大向量尺寸 returns 错误的结果

Question

我正在尝试实现一个 CUDA 内核，它正在计算两个向量的点积。对于小向量大小，代码工作正常，我得到正确的结果，但对于更大的向量，它只是计算错误。我正在实施三种不同的方法来计算点积：

serial版本（直接在c++中，没有任何优化）
CUDA 内核
CUBLAS版本

我在 cpp 文件中的主要内容如下所示：

float *h_x,*h_y;
float res1=0.0, res2=0.0, res3=0.0;

h_x=(float*)malloc(Main::N*sizeof(float)); random_ints_Vec(h_x);
h_y=(float*)malloc(Main::N*sizeof(float)); random_ints_Vec(h_y);

double serialTimer;
double cublasTimer;
double cudaTimer;

res1=serial_dotProd(h_x,h_y,&serialTimer);  
res2=cublas_dotProd(h_x,h_y,&cublasTimer);
res3=cuda_dotProd(h_x,h_y,&cudaTimer);      

free(h_x); free(h_y);

序列号版本：

float Main::serial_dotProd(float* x, float* y, double* time){
std::clock_t start;
start = std::clock();

float res1=0.0;
for (int i=0;i<Main::N;++i) {
    res1+=x[i]*y[i];
}

*time= ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
return res1;}

CUDA版本：

float Main::cuda_dotProd(float *h_x,float *h_y,double* time){
float *d_x,*d_y,*d_res, *h_res;
h_res=(float*)malloc(Main::BLOCKS_PER_GRID*sizeof(float));

size_t bfree, afree, total;
cudaMemGetInfo(&bfree,&total);

cudaMalloc((void**) &d_x,Main::N*sizeof(float));
cudaMalloc((void**) &d_y,Main::N*sizeof(float));
cudaMalloc((void**) &d_res,Main::BLOCKS_PER_GRID*sizeof(float));
cudaCheckErrors("cuda malloc fail");

cudaMemGetInfo(&afree,&total);
std::cout<<" > memory used for cuda-version:"<< (bfree -afree)/1048576<< "MB ("<<total/1048576 <<"MB)" <<"\n";

cudaMemcpy(d_x,h_x,Main::N*sizeof(float),cudaMemcpyHostToDevice);
cudaMemcpy(d_y,h_y,Main::N*sizeof(float),cudaMemcpyHostToDevice);   

std::clock_t start;
start = std::clock();   
DotProdWrapper(d_x,d_y,d_res,(Main::N+Main::THREADS_PER_BLOCK-1)/Main::THREADS_PER_BLOCK,Main::THREADS_PER_BLOCK,Main::N);

*time= ( std::clock() - start ) / (double) CLOCKS_PER_SEC;

cudaMemcpy(h_res,d_res,Main::BLOCKS_PER_GRID*sizeof(float),cudaMemcpyDeviceToHost);

float c= 0;
for (int i=0;i<Main::BLOCKS_PER_GRID;++i){
  c+=h_res[i];
}
cudaFree(d_x);
cudaFree(d_y);
cudaFree(d_res);    

free(h_res);
return c;}

CUDA 内核：

__global__ void DotProd(float* x, float* y, float* scalar,unsigned long int N){
    extern __shared__ float cache[];

    int tid = threadIdx.x+ blockIdx.x*blockDim.x;
    int cacheIndex = threadIdx.x;

    float temp=0;
    while (tid<N){
        temp+=x[tid]*y[tid];
        tid +=blockDim.x*gridDim.x; 
    }
    cache[cacheIndex]=temp;
    __syncthreads();

    int i=blockDim.x/2;
    while(i!=0){
        if (cacheIndex<i)
            cache[cacheIndex]+=cache[cacheIndex+i];
        __syncthreads();
        i/=2;
    }
    if(cacheIndex==0)
        scalar[blockIdx.x]=cache[cacheIndex];
}

CUBLAS 版本：

float Main::cublas_dotProd(float *h_x,float *h_y, double* time){
    float *d_x,*d_y;
    float *res;
    float result=0.0;
    cublasHandle_t h;
    cublasCreate(&h);
    cublasSetPointerMode(h, CUBLAS_POINTER_MODE_DEVICE);

    size_t bfree, afree, total;
    cudaMemGetInfo(&bfree,&total);

    cudaMalloc((void**) &d_x,Main::N*sizeof(float));
    cudaMalloc((void**) &d_y,Main::N*sizeof(float));
    cudaMalloc( (void **)(&res), sizeof(float) );
    cudaCheckErrors("cublas malloc fail");

    cudaMemGetInfo(&afree,&total);

     cudaMemcpy(d_x, h_x, Main::N*sizeof(float), cudaMemcpyHostToDevice);
     cudaMemcpy(d_y, h_y, Main::N*sizeof(float), cudaMemcpyHostToDevice);

    cublasSetVector(Main::N,sizeof(float),h_x,1,d_x,1);
    cublasSetVector(Main::N,sizeof(float),h_y,1,d_y,1);


    std::clock_t start;
    start = std::clock();

    cublasSdot(h,Main::N,d_x,1,d_y,1,res);

    *time= ( std::clock() - start ) / (double) CLOCKS_PER_SEC;

    cudaMemcpy(&result, res, sizeof(float), cudaMemcpyDeviceToHost);

    cudaFree(d_x);
    cudaFree(d_y);
    cudaFree(res);
    return result;
}

我在不同设置下计算得到的结果如下：

N=204800 , THREADS_PER_BLOCK:512, BLOCKS_PER_GRID:400 serial_dotProd=4.15369e+06 ; cublas_dotProd=4.15369e+06 ; cuda_dotProd=4.15369e+06
N=20480000 , THREADS_PER_BLOCK:512, BLOCKS_PER_GRID:40000 serial_dotProd=4.04149e+08 ; cublas_dotProd=4.14834e+08 ; cuda_dotProd=4.14833e+08

我不知道为什么，但在我的向量达到一定大小后，我得到了错误的结果。矢量确实适合 SDRAM，每个块的共享内存也有足够的 space 来分配内存。提前致谢。

Answer 1

这个问题经常出现，以至于 Nvidia 专门 an entire section of the CUDA Floating Point and IEEE 754 guide to it. It's also briefly mentioned in the CUDA C Best Practices Guide。

简短的解释是双重的：

与相应的精确数学运算不同，由于涉及舍入误差，浮点算术运算不具有关联性。这意味着将求和从直接串行求和重新排序为适合并行执行的树结构将稍微改变结果（随着求和的值数量的增加更是如此）。巧合的是，在大多数情况下，树状排列也给出了比顺序求和更接近精确数学求和的结果。
CUDA 编译器倾向于更积极地使用融合乘加（FMA，一种先乘后加运算，其中省略了中间舍入阶段）。同样，数学上正确的结果往往更接近使用 FMA 获得的结果。

所以可能的答案是使用 CUDA 获得的结果可能比简单的串行 CPU 版本更接近准确结果（这就是为什么我要求您使用更高的精度再次执行实验）。

带有 CUDA 内核的点积用于大向量尺寸 returns 错误的结果

Dot Product with a CUDA kernel for big vector sizes returns wrong results

c++

cuda

cublas