非方阵的二维内核调用和启动参数

Posted 2023-03-23

【中文标题】非方阵的二维内核调用和启动参数

【英文标题】：2D kernel calling and launch parameters for non-square matrix

【发布时间】：2012-08-02 14:18:08

【问题描述】：

我正在尝试将以下（简化的）嵌套循环移植为 CUDA 2D 内核。 NgS 和 NgO 的大小会随着更大的数据集而增加；现在我只想让这个内核为所有值输出正确的结果：

// macro that translates 2D [i][j] array indices to 1D flattened array indices
#define idx(i,j,lda) ( (j) + ((i)*(lda)) )

int NgS  = 1859;
int NgO  = 900;

// 1D flattened matrices have been initialized as:
 Radio_cpu = new double [NgS*NgO];
Result_cpu = new double [NgS*NgO];
// ignoring the part where they are filled w/ data

for (m=0; m<NgO; m++)         
    for (n=0; n<NgS; n++) 
            Result_cpu[idx(n,m,NgO)]] = k0*Radio_cpu[idx(n,m,NgO)]];
    

我遇到的示例通常处理方形循环，与 CPU 版本相比，我无法获得所有 GPU 数组索引的正确输出。这是调用内核的主机代码：

dim3 dimBlock(16, 16);
dim3 dimGrid;
dimGrid.x = (NgO + dimBlock.x - 1) / dimBlock.x;
dimGrid.y = (NgS + dimBlock.y - 1) / dimBlock.y;

// Result_gpu and Radio_gpu are allocated versions of the CPU variables on GPU
trans<<<dimGrid,dimBlock>>>(NgO, NgS, k0, Radio_gpu, Result_gpu);

这是内核：

__global__ void trans(int NgO, int NgS,
                      double k0, double * Radio, double * Result) 

int n = blockIdx.x * blockDim.x + threadIdx.x;
int m = blockIdx.y * blockDim.y + threadIdx.y;

if(n > NgS || m > NgO) return;

// map the two 2D indices to a single linear, 1D index
int grid_width = gridDim.x * blockDim.x;
int idxxx = m + (n * grid_width);

Result[idxxx] = k0 * Radio[idxxx];

使用当前代码，我继续比较 Result_cpu 变量和 Result_gpu 变量复制回来。当我循环浏览我得到的值时：

    // matches from NgS = 0...913
    Result_gpu[NgS = 913][NgO = 0]: -56887.2
    Result_cpu[Ngs = 913][NgO = 0]: -56887.2

    // mismatches from NgS = 914...1858
    Result_gpu[NgS = 914][NgO = 0]: -12.2352
    Result_cpu[NgS = 914][NgO = 0]: 79448.6

无论NgO 的值如何，此模式都是相同的。我一直试图通过查看各种示例并尝试更改来找出我犯了错误的地方，但到目前为止，这个方案已经有效地减去了手头的明显问题，而其他方案已经导致内核调用错误/离开GPU 数组未对所有值进行初始化。由于我显然看不到错误，如果有人能指出我正确的方向来解决问题，我将不胜感激。我很确定它就在我的鼻子下面，我看不到它。

以防万一，我正在 Kepler 卡上测试此代码，使用 MSVC 2010、CUDA 4.2 和 304.79 驱动程序进行编译，并使用 arch=compute_20,code=sm_20 和 arch=compute_30,code=compute_30 标志编译代码，没有区别。

【参考方案1】：

@vaca_loca：我测试了以下内核（它也适用于非方形块尺寸）：

__global__ void trans(int NgO, int NgS,
                  double k0, double * Radio, double * Result) 

int n = blockIdx.x * blockDim.x + threadIdx.x;
int m = blockIdx.y * blockDim.y + threadIdx.y;
if(n > NgO || m > NgS) return;
int ofs = m * NgO + n;
Result[ofs] = k0 * Radio[ofs];


void test() 

int NgS  = 1859, NgO  = 900;
int data_sz = NgS * NgO, bytes = data_sz * sizeof(double);
cudaSetDevice(0);
double *Radio_cpu = new double [data_sz*3],
    *Result_cpu = Radio_cpu + data_sz,
    *Result_gpu = Result_cpu + data_sz;
double k0 = -1.7961233;

srand48(time(NULL));
int i, j, n, m;
for(m=0; m<NgO; m++) 
  for (n=0; n<NgS; n++) 
        Radio_cpu[m + n*NgO] = lrand48() % 234234;
        Result_cpu[m + n*NgO] = k0*Radio_cpu[m + n*NgO];
    


double *g_Radio, *g_Result;
cudaMalloc((void **)&g_Radio, bytes * 2);
g_Result = g_Radio + data_sz;
cudaMemcpy(g_Radio, Radio_cpu, bytes, cudaMemcpyHostToDevice);

dim3 dimBlock(16, 16);
dim3 dimGrid;
dimGrid.x = (NgO + dimBlock.x - 1) / dimBlock.x;
dimGrid.y = (NgS + dimBlock.y - 1) / dimBlock.y;

trans<<<dimGrid,dimBlock>>>(NgO, NgS, k0, g_Radio, g_Result);

cudaMemcpy(Result_gpu, g_Result, bytes, cudaMemcpyDeviceToHost);

for(m=0; m<NgO; m++) 
    for (n=0; n<NgS; n++) 
        double c1 = Result_cpu[m + n*NgO],
                c2 = Result_gpu[m + n*NgO];
        if(std::abs(c1-c2) > 1e-4)
            printf("(%d;%d): %.7f %.7f\n", n, m, c1, c2);
    

cudaFree(g_Radio);
delete []Radio_cpu;

不过，在我看来，使用四边形访问全局内存中的数据可能对缓存不太友好，因为访问步幅非常大。如果您的算法在 2D 局部性中访问数据至关重要，您可能会考虑使用 2D 纹理

【讨论】：

我还以为是这样，但是当我在条件语句中交换变量时，我收到内核启动错误，因此所有 Result_GPU 值都是 0。

我正在尝试使用内核中的 printf 语句进行调试，以查看同时索引如何出错。我觉得可能存在多个问题——我可能没有正确调用内核并且我可能没有正确地在内核中建立索引。

我编辑了我的答案：再看一遍，这对我来说很好

您的 ofs 定义是我尝试的第一件事，但它对我不起作用——当我添加此示例中的更改时，它开始主要起作用：code.google.com/p/stanford-cs193g-sp2010/wiki/…

我想我可能已经通过运行该示例几次来解决问题 - 似乎 blocksize.x 和 blocksize.y 应该与相应的 num_elements_x 和 num_elements_y 有一个共同因素，但我这个想法还是要在实际代码中测试一下。