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

问题描述

我尝试实现某种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,cudaMemcpyHostToDevice);
    cudaMemcpy(d_phi1,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;
}

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

解决方法

在CUDA中，当指定kernel execution configuration arguments时，首先是网格配置，然后是块配置。

因此，给定变量dimGrid和dimBlock的定义和用法，这是不正确的：

    Jacobi<<<dimBlock,dimGrid>>>

应该是：

    Jacobi<<<dimGrid,dimBlock>>>

网格和块配置变量具有不同的限制，因此，由于混合违反了硬件限制，在520,520的问题维度上发生了第一次启动失败