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

Question

在我进入之前，这是关于正在发生的事情的一般想法：

总的想法是，我有x个浮点数数组，我希望将每个数组依次添加到另一个数组中(标量添加)：

T=数组；

A=数组；

T=零

T += a

T += a1

..。

T += aN

其中+=表示标量加法。

这是直接的。我试图缩小代码，我必须尽可能紧凑，并保留功能。这里的问题是，对于某些大小的数组，我看到任何大于128x128x108的问题。基本上，复制回主机的内存的总和与我计算出来的不一样。我被困在这上面一整天了，所以我不会再浪费时间了。我真的无法解释为什么会发生这种事。我的理由是：

• 使用太多的常量空间(不使用任何)
• 使用过多的寄存器(否)
• 内核中检查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: %d\n",e);
        printf("%s\n",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 depth\n");
        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 : %d\n",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 allocated\n");

    // 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 value\n");
                exit(0);
            }
            count++;
        }
        printf("Total Count: %d\n",count);
        printf("Real Array sum: %18f\nTotal values counted : %d\n",tempsum,count*(i+1));
        printf("Calculated Array sum: %ld\n\n",(i+1)*fMemSize );
    }

    for(int i=0; i<nHostF; i++){
        cudaFree(f_dev[i]);
    }

    cudaFree(f_temp_dev);
    printf("Memory free. Program successfully complete\n");
    delete f_dev;
    delete f_temp_host;
}

Answer 1

你的设备代码没有什么问题。正在发生的一切是，在较大的问题大小时，您正在耗尽单精度浮点的能力，以精确计算代码在大运行大小下生成的大整数值。如果将主机端求和代码替换为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 value\n");
            exit(0);
        }
        count++;
    }

您应该发现代码在更大的大小下按预期运行。或者，可以用双精度算法来完成主机端的求和。如果你还没有读过，我强烈推荐http://download.oracle.com/docs/cd/E19957-01/806-3568/ncg_goldberg.html。它将有助于解释您在这个示例中出错的地方，并且它所赋予的智慧可能有助于防止将来发生类似的错误。