c++ - 大数组的C++/CUDA怪异行为

Question

我尝试实现某种Jacobi算法，并测量不同网格大小所花费的时间。
对于具有相同数量的迭代，无论网格有多大，我都不使用某种形式的残差，而是让算法始终运行4000个迭代(但数组的大小不同)。
直到我超过510x510网格(两倍)为止，它完全可以正常工作。 510x510大约需要2763498微秒，然后520x520大约需要1778微秒。
我已经尝试从双数组更改为浮点数组，以确保这不是某种形式的内存不足，但我不知道问题出在哪里。

__global__ void Jacobi(double *a, double *b, double *c, int L){
int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;

if(row > 0 && col > 0 && row < L-1 && col < L-1){
    a[row * L + col] = 1.0/4.0 * (b[col+1 + row*L] + b[col - 1 + row*L] + b[col + (row+1)*L] + b[col + (row-1)*L] - c[col + row*L]);
    __syncthreads();
    b[row*L + col] = a[row*L+col];
    }
}

int main(){
int L;
int Iterations;
double *h_phi1;
double *h_phi2;

double *h_f;
FILE * temp = fopen("Timings.out", "w");
for (L=10;L<10000;L+=10){
    long long int size = L*L*sizeof(double);
    h_f = (double*) malloc(size);
    h_phi1 = (double*) malloc(size);
    h_phi2 = (double*) malloc(size);


    for(int i=0;i<L;i++){
        for(int j=0;j<L;j++){
            h_f[j+i*L] = (pow(1.0/float(L),2.0))*exp(-100.0*(pow((float(i)/float(L) - float(1.0/3.0) ),2.0) + pow(( float(j)/float(L) - float(1.0/3.0) ),2.0))) - 
            (pow(1.0/ float(L),2.0))*exp(- 100.0*(pow(( float(i)/ float(L) -(1.0- 1.0/3.0)),2.0) + pow(( float(j)/float(L)-(1.0-float(1.0/3.0))),2.0)));
            h_phi1[j+i*L] = 0;
            h_phi2[j+i*L] = 0;
        }
    }

    //allocate memory on GPU
    double *d_phi1;
    cudaMalloc(&d_phi1, size);
    double *d_phi2;
    cudaMalloc(&d_phi2, size);
    double *d_f;
    cudaMalloc(&d_f, size);
    
    //set CTA
    int threads = 16;
    int blocks = (L+threads-1)/threads;
    double epsCalc;

    //Setup Kernel launch parameters
    dim3 dimBlock(threads, threads);
    dim3 dimGrid(blocks, blocks);

    //Setup timing and Cpy Memory from Host to Device
    Iterations = 0;
    auto t1 = std::chrono::high_resolution_clock::now();
    cudaMemcpy(d_phi2, h_phi2, size, cudaMemcpyHostToDevice);
    cudaMemcpy(d_f, h_f, size, cudaMemcpyHostToDevice);
    cudaMemcpy(d_phi1, h_phi2, size, cudaMemcpyHostToDevice);

    //Launch Kernel
    for (int j=0;j<4000;j++){
        Iterations += 1;
        Jacobi<<<dimBlock, dimGrid>>>(d_phi2,d_phi1,d_f,L);

    }
    auto t2 = std::chrono::high_resolution_clock::now();
    auto duration = std::chrono::duration_cast<std::chrono::microseconds>( t2 - t1 ).count();
    
    fprintf(temp, "%lf % %d \n", L, duration);
    printf("I reached the end of Jacobi after %d Iterations!\n Time taken was %d in milliseconds", Iterations, duration);
    cudaFree(d_f); cudaFree(d_phi2), cudaFree(d_phi1);
    free(h_f); free(h_phi2); free(h_phi1);
}

return 0;
}

我希望有人可以引导我到犯错的地方。

Answer 1

在CUDA中，当指定kernel execution configuration arguments时，首先是网格配置，然后是块配置。
因此，给定变量dimGrid和dimBlock的定义和用法，这是不正确的:

    Jacobi<<<dimBlock, dimGrid>>>

它应该是:

    Jacobi<<<dimGrid, dimBlock>>>

网格和块配置变量具有不同的限制，因此，由于混合违反了硬件限制(发生在问题尺寸为520,520)而导致的第一次启动失败发生