大数组的 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

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

C++/CUDA weird behavior with large arrays

c++

algorithm

cuda

large-data