Cuda内核中的大型for循环不适用于大型数组[关闭]

Question

我已经使用 Cuda 实现了各种算法，例如矩阵乘法、Cholesky 分解和下三角矩阵的求逆（通过前向替换）。

对于其中一些算法，我在内核中有一个 for 循环，它会多次重复部分内核代码。这一切都适用于（扁平化：由一维数组表示）矩阵（浮点数）高达约 200x200，for 循环调用内核代码的一部分 200 次。将矩阵大小增加到 1000x1000（使用 for 循环调用内核代码的部分 1000 次）会使 GPU 花费尽可能多的计算时间，这取决于使用较小矩阵大小的试验。但是似乎没有运行内核代码（包括 for 循环之外的部分）（输出矩阵自初始化以来没有任何元素发生变化）。如果我将矩阵大小增加到 500 左右，如果我将 for 循环中的限制器设置为某个低值（例如 3），我有时可以让内核运行。

我是否在这里遇到了一些硬件限制，或者有什么技巧可以让这些 for 循环适用于大型矩阵？

这是一个完整的代码示例，您可以将其复制到 .cu 文件中。内核尝试将矩阵 A (W*H) 的内容复制到矩阵 B (W*H)。输出显示了两个矩阵的第一个元素，对于 W*H

  #include <cuda.h>
  #include <cuda_runtime.h>
  #include <cuda_runtime_api.h>
  #include <stdio.h> 
  #include <stdlib.h>
  #include <math.h>
  #include <iostream>
  #include <time.h>

#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
inline void gpuAssert(cudaError_t code, char *file, int line, bool abort=true)
{
   if (code != cudaSuccess) 
   {
      fprintf(stderr,"GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
      if (abort) exit(code);
   }
}  

__global__ void MatrixCopy(float *A, float *B, int W)
{

int i = blockIdx.x*blockDim.x + threadIdx.x;
int j = blockIdx.y*blockDim.y + threadIdx.y;

B[j*W + i]=A[j*W + i];

}

int main(void)
{

clock_t start1=clock();

int W=1000;
int H=1000;
float *A, *B;
float *devA, *devB;

A=(float*)malloc(W*H*sizeof(float));
B=(float*)malloc(W*H*sizeof(float));

for(int i=0; i<=W*H; i++)
{
    A[i]=rand() % 3;
    A[i]=A[i]+1;
    B[i]=0;
}

gpuErrchk( cudaMalloc( (void**)&devA, W*H*sizeof(float) ) ); 
gpuErrchk( cudaMalloc( (void**)&devB, W*H*sizeof(float) ) ); 

gpuErrchk( cudaMemcpy( devA, A, W*H*sizeof(float), cudaMemcpyHostToDevice ) );
gpuErrchk( cudaMemcpy( devB, B, W*H*sizeof(float), cudaMemcpyHostToDevice ) );

dim3 threads(32,32);
int bloW=(int)ceil((double)W/32);
int bloH=(int)ceil((double)H/32);
dim3 blocks(bloW, bloH);

clock_t finish1=clock();
clock_t start2=clock();

MatrixCopy<<<blocks,threads>>>(devA, devB, W);
gpuErrchk( cudaPeekAtLastError() );

gpuErrchk( cudaMemcpy( B, devB, W*H*sizeof(float), cudaMemcpyDeviceToHost ) );

clock_t finish2=clock();

printf("\nGPU calculation time (ms): %d\nInitialization time (ms): %d\n\n", (int)ceil(double(((finish2-start2)*1000/(CLOCKS_PER_SEC)))), (int)ceil(double(((finish1-start1)*1000/(CLOCKS_PER_SEC)))));
printf("\n%f\n", A[0]);
printf("\n%f\n\n", B[0]);

gpuErrchk( cudaFree(devA) );
gpuErrchk( cudaFree(devB) );

free(A);
free(B);

#ifdef _WIN32 
    system ("PAUSE"); 
#endif 

return 0;

}

Answer 1

您的内核没有线程检查。

您正在像这样决定网格大小（以块为单位）：

int bloW=(int)ceil((double)W/32);
int bloH=(int)ceil((double)H/32);

对于H 和W 的值，它们甚至不是每个块大小的线程数的倍数 (32)，这会在您关心的实际矩阵 (1000x1000) 之外创建额外的线程和块。这没有什么问题；这是常见的做法。

但是，我们必须确保那些额外的线程实际上不做任何事情（即不产生对内存的无效访问）。您的内核不提供此检查。

如果你修改你的内核是这样的：

__global__ void MatrixCopy(float *A, float *B, int W, int H)
{

  int i = blockIdx.x*blockDim.x + threadIdx.x;
  int j = blockIdx.y*blockDim.y + threadIdx.y;

  if ((i <  W) && (j < H))
    B[j*W + i]=A[j*W + i];

}

我认为你会有更好的结果。没有这个，你在内核中的一些A 和B 引用会产生越界访问，如果你用cuda-memcheck 运行你的代码，你可以看到。您还必须修改内核调用行以添加H 参数。我还没有真正弄清楚你的i变量是对应H还是W；我假设您可以这样做并在需要时进行更改。在这种情况下，由于矩阵是方形的，所以并不重要。

您应该在遇到 CUDA 代码问题时使用proper cuda error checking。我建议您在此处发帖寻求帮助之前执行此操作。