返回错误值的数组的顺序求和

Question

在我进入之前，这是对正在发生的事情的总体思路：

一般的想法是我有 x 个浮点数组，我想将每个数组按顺序添加到另一个数组（标量添加）：

t = 数组；

a = 数组数组；

t = 零

t += a[0]

t += a[1]

...

t += a[N]

其中 += 表示标量加法。

这是直截了当的。我试图缩小我必须尽可能紧凑并保留功能的代码。这里的问题是，对于某些大小的数组 - 我看到任何大于 128 x 128 x 108 的问题。基本上复制回主机的内存总和与我计算的不一样。我整天都被困在这个问题上，所以我要停止浪费我的时间了。我真的无法解释为什么会这样。我已经推理了：

• 使用过多的常量空间（不使用任何空间）
• 使用过多的寄存器（否）
• 内核中检查 idx、idy、idz 是否在界限内的条件不正确（可能仍然如此）
• gpu 有点有趣（在 gt280 以及 tesla C1060 和 C2060 上试过）
• 不正确的 printf 格式（希望是这样） *...

该列表可以继续。如果您有时间，感谢您浏览此内容。该问题似乎与内存有关（即 > 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 summation替换你的主机端求和代码，像这样：

    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++;
    }

您应该会发现代码在较大尺寸下运行正常。或者，主机端求和可以用双精度算术代替。如果你还没有读过，我强烈推荐What Every Computer Scientist Should Know About Floating-Point Arithmetic。它将有助于解释您在此示例中出错的地方，并且它所传授的智慧可能有助于防止将来犯下类似的faux pas。