返回不正确值的数组的顺序求和

sequential summation of arrays returning incorrect value

本文关键字:顺序 求和 数组 返回 不正确      更新时间:2023-10-16

在我进入之前,这是正在发生的事情的总体思路:

一般的想法是,我有x个数组的浮点数,我想把每个数组依次添加到另一个数组(标量添加):

a = array;

t = 0

t += a[0]

t += a[1]

t += a[N]

其中+=表示标量加法。

这很直接。我试图缩减代码,我必须尽可能紧凑,并保留功能。这里的问题是,对于某些大小的数组-我看到任何大于128 x 128 x 108的问题。基本上,复制回主机的内存总和和我计算出来的不一样。我一整天都在纠结这件事,所以我不想再浪费时间了。我真的无法解释为什么会这样。

  • 使用太多常量空间(不使用任何)
  • 使用太多寄存器(no)
  • 内核检查idx,idy,idz是否在边界内的错误条件(这仍然可能是它)
  • 一些有趣的gpu(尝试在gt280,和特斯拉C1060和C2060)
  • 打印格式不正确(我希望是这样)*…

这样的例子不胜枚举。如果你有时间,谢谢你浏览这个。这个问题似乎与内存有关(即> 128*128*108的内存大小不起作用)。所以64*128*256 可以工作,或者它们的任何排列)。

这是完整的源代码(应该可以用nvcc编译):

#include <cuda.h>
#include <iostream>
#include <stdio.h>
#include <assert.h>
#define BSIZE 8
void cudaCheckError(cudaError_t e,const char * msg) {
    if (e != cudaSuccess){
        printf("Error number: %dn",e);
        printf("%sn",msg);
    }
};
__global__ void accumulate(float * in,float * out, int3 gdims, int zlevel) {
    int idx = blockIdx.x*blockDim.x + threadIdx.x;
    int idy = blockIdx.y*blockDim.y + threadIdx.y;
    int idz = threadIdx.z;
    long int index = (zlevel*((int)BSIZE)+idz)*gdims.x*gdims.y+ 
        idy*gdims.x+ 
        idx;
    if ( idx < gdims.x && idy < gdims.y && (idz + zlevel*(int)BSIZE) < gdims.z) {
        out[index] += in[index];
    }
};
int main(int argc, char * argv[]) {
    int width, 
    height,
    depth; 
    if (argc != 4) {
        printf("Must have 3 inputs: width height depthn");
        exit(0);
    }
    float tempsum;
    int count =0;
    width = atoi(argv[1]);
    height = atoi(argv[2]);
    depth = atoi(argv[3]);
    printf("Dimensions (%d,%d,%d)n",width,height,depth);
    int3 dFull;
    dFull.x = width+2;
    dFull.y = height+2;
    dFull.z = depth+2;
    printf("Dimensions (%d,%d,%d)n",dFull.x,dFull.y,dFull.z);
    int fMemSize=dFull.x*dFull.y*dFull.z;
    int nHostF=9;
    float * f_hostZero;
    float ** f_dev;
    float * f_temp_host;
    float * f_temp_dev;
    dim3 grid( dFull.x/(int)BSIZE+1, dFull.y/(int)BSIZE + 1);
    dim3 threads((int)BSIZE,(int)BSIZE,(int)BSIZE);
    printf("Threads (x,y) : (%d,%d)nGrid (x,y) : (%d,%d)n",threads.x,threads.y,grid.x,grid.y);
    int num_zsteps=dFull.z/(int)BSIZE + 1;
    printf("Number of z steps to take : %dn",num_zsteps);
    // Host array allocation
    f_temp_host = new float[fMemSize];
    f_hostZero = new float[fMemSize];

    // Allocate nHostF address on host 
    f_dev = new float*[nHostF];
    // Host array assignment
    for(int i=0; i < fMemSize; i++){
        f_temp_host[i] = 1.0;
        f_hostZero[i] = 0.0;
    }
    // Device allocations - allocated for array size + 2
    for(int i=0; i<nHostF; i++){
        cudaMalloc((void**)&f_dev[i],sizeof(float)*fMemSize);
    }

    // Allocate the decive pointer
    cudaMalloc( (void**)&f_temp_dev, sizeof(float)*fMemSize);
    cudaCheckError(cudaMemcpy((void *)f_temp_dev,(const void *)f_hostZero,
        sizeof(float)*fMemSize,cudaMemcpyHostToDevice),"At first mem copy");
    printf("Memory regions allocatedn");
    // Copy memory to each array
    for(int i=0; i<nHostF; i++){
        cudaCheckError(cudaMemcpy((void *)(f_dev[i]),(const void *)f_temp_host,
            sizeof(float)*fMemSize,cudaMemcpyHostToDevice),"At first mem copy");
    }
    // Add value 1.0 (from each array n f_dev[i]) to f_temp_dev
    for (int i=0; i<nHostF; i++){
        for (int zLevel=0; zLevel<num_zsteps; zLevel++){
            accumulate<<<grid,threads>>>(f_dev[i],f_temp_dev,dFull,zLevel);
            cudaThreadSynchronize();
        }
        cudaCheckError(cudaMemcpy((void *)f_temp_host,(const void *)f_temp_dev,
            sizeof(float)*fMemSize,cudaMemcpyDeviceToHost),"At mem copy back");
        tempsum=0.f;
        count =0;
        for(int k = 0 ; k< fMemSize; k++){
            tempsum += f_temp_host[k];
            assert ( (int)f_temp_host[k] == (i+1) );
            if ( f_temp_host[k] !=(float)(i+1) ) {
                printf("Found invalid return valuen");
                exit(0);
            }
            count++;
        }
        printf("Total Count: %dn",count);
        printf("Real Array sum: %18fnTotal values counted : %dn",tempsum,count*(i+1));
        printf("Calculated Array sum: %ldnn",(i+1)*fMemSize );
    }
    for(int i=0; i<nHostF; i++){
        cudaFree(f_dev[i]);
    }
    cudaFree(f_temp_dev);
    printf("Memory free. Program successfully completen");
    delete f_dev;
    delete f_temp_host;
}

您的设备代码没有问题。所发生的一切是,在大问题规模时,您正在耗尽单精度浮点数的能力,以精确计算代码在大运行规模下产生的大整数值。如果将主机端求和代码替换为Kahan求和代码,如下所示:

    tempsum=0.f;
    count =0;
    float c=0.f;
    for(int k = 0 ; k< fMemSize; k++){
        float y = f_temp_host[k] - c;
        float t = tempsum + y;
        c = (t - tempsum) - y;
        tempsum = t;
        assert ( (int)f_temp_host[k] == (i+1) );
        if ( f_temp_host[k] !=(float)(i+1) ) {
            printf("Found invalid return valuen");
            exit(0);
        }
        count++;
    }

时,您应该会发现代码在更大的尺寸下按预期运行。或者,主机端求和可以用双精度算术来代替。如果你还没有读过它,我强烈推荐每个计算机科学家应该知道的浮点算术。它将帮助你解释在这个例子中哪里出错了,它所赋予的智慧可能有助于防止在将来犯类似的错误