cuda - CUDA 中的二维时域有限差分 (FDTD)

Question

我想用 CUDA 语言编写电磁二维时域有限差分 (FDTD) 代码。 磁场更新的C代码如下

// --- Update for Hy and Hx
for(int i=n1; i<=n2; i++)
   for(int j=n11; j<=n21; j++){
      Hy[i*ydim+j]=A[i*ydim+j]*Hy[i*ydim+j]+B[i*ydim+j]*(Ezx[(i+1)*ydim+j]-Ezx[i*ydim+j]+Ezy[(i+1)*ydim+j]-Ezy[i*ydim+j]);
  Hx[i*ydim+j]=G[i*ydim+j]*Hx[i*ydim+j]-H[i*ydim+j]*(Ezx[i*ydim+j+1]-Ezx[i*ydim+j]+Ezy[i*ydim+j+1]-Ezy[i*ydim+j]);
   }
}

我的第一个并行化尝试是以下内核:

__global__ void H_update_kernel(double* Hx_h, double* Hy_h, double* Ezx_h, double* Ezy_h, double* A_h, double* B_h,double* G_h, double* H_h, int n1, int n2, int n11, int n21)
{
   int idx = blockIdx.x*BLOCK_SIZE_X + threadIdx.x;
   int idy = blockIdx.y*BLOCK_SIZE_Y + threadIdx.y;

   if ((idx <= n2 && idx >= n1)&&(idy <= n21 && idy >= n11)) {
      Hy_h[idx*ydim+idy]=A_h[idx*ydim+idy]*Hy_h[idx*ydim+idy]+B_h[idx*ydim+idy]*(Ezx_h[(idx+1)*ydim+idy]-Ezx_h[idx*ydim+idy]+Ezy_h[(idx+1)*ydim+idy]-Ezy_h[idx*ydim+idy]);
  Hx_h[idx*ydim+idy]=G_h[idx*ydim+idy]*Hx_h[idx*ydim+idy]-H_h[idx*ydim+idy]*(Ezx_h[idx*ydim+idy+1]-Ezx_h[idx*ydim+idy]+Ezy_h[idx*ydim+idy+1]-Ezy_h[idx*ydim+idy]); }

}

但是，通过同时使用 Visual Profiler，我对这个解决方案并不满意，原因有二: 1)内存访问合并不佳； 2) 没有使用共享内存。

然后我决定使用以下解决方案

__global__ void H_update_kernel(double* Hx_h, double* Hy_h, double* Ezx_h, double* Ezy_h, double* A_h, double* B_h,double* G_h, double* H_h, int n1, int n2, int n11, int n21)
{
    int i       = threadIdx.x;
int j       = threadIdx.y;
int idx     = blockIdx.x*BLOCK_SIZE_X + threadIdx.x;
int idy     = blockIdx.y*BLOCK_SIZE_Y + threadIdx.y;

int index1  = j*BLOCK_SIZE_Y+i;

int i1      = (index1)%(BLOCK_SIZE_X+1);
int j1      = (index1)/(BLOCK_SIZE_Y+1);

int i2      = (BLOCK_SIZE_X*BLOCK_SIZE_Y+index1)%(BLOCK_SIZE_X+1);
int j2      = (BLOCK_SIZE_X*BLOCK_SIZE_Y+index1)/(BLOCK_SIZE_Y+1);

__shared__ double Ezx_h_shared[BLOCK_SIZE_X+1][BLOCK_SIZE_Y+1];     
__shared__ double Ezy_h_shared[BLOCK_SIZE_X+1][BLOCK_SIZE_Y+1];     

if (((blockIdx.x*BLOCK_SIZE_X+i1)<xdim)&&((blockIdx.y*BLOCK_SIZE_Y+j1)<ydim))
    Ezx_h_shared[i1][j1]=Ezx_h[(blockIdx.x*BLOCK_SIZE_X+i1)*ydim+(blockIdx.y*BLOCK_SIZE_Y+j1)];

if (((i2<(BLOCK_SIZE_X+1))&&(j2<(BLOCK_SIZE_Y+1)))&&(((blockIdx.x*BLOCK_SIZE_X+i2)<xdim)&&((blockIdx.y*BLOCK_SIZE_Y+j2)<ydim)))
    Ezx_h_shared[i2][j2]=Ezx_h[(blockIdx.x*BLOCK_SIZE_X+i2)*xdim+(blockIdx.y*BLOCK_SIZE_Y+j2)];

__syncthreads();

if ((idx <= n2 && idx >= n1)&&(idy <= n21 && idy >= n11)) {
    Hy_h[idx*ydim+idy]=A_h[idx*ydim+idy]*Hy_h[idx*ydim+idy]+B_h[idx*ydim+idy]*(Ezx_h_shared[i+1][j]-Ezx_h_shared[i][j]+Ezy_h[(idx+1)*ydim+idy]-Ezy_h[idx*ydim+idy]);
    Hx_h[idx*ydim+idy]=G_h[idx*ydim+idy]*Hx_h[idx*ydim+idy]-H_h[idx*ydim+idy]*(Ezx_h_shared[i][j+1]-Ezx_h_shared[i][j]+Ezy_h[idx*ydim+idy+1]-Ezy_h[idx*ydim+idy]); }

    } 

需要使用索引技巧来使 BS_x * BS_y 线程 block 将 (BS_x+1)*(BS_y+1) 个全局内存位置读取到共享内存。 我相信这个选择比前一个更好，因为使用了共享内存，虽然并不是所有的访问都真正合并，请参阅

Analyzing memory access coalescing of my CUDA kernel

我的问题是，如果你们中的任何人都可以向我提出在合并内存访问方面更好的解决方案。谢谢。

Answer 1

感谢您提供分析信息。

您的第二个算法更好，因为您获得了更高的 IPC。尽管如此，在 CC 2.0 上，最大 IPC 为 2.0，因此您在第二个解决方案中的平均值为 1.018 意味着仅使用了一半的可用计算能力。通常，这意味着您的算法是内存限制的，但我不确定您的情况，因为内核中的几乎所有代码都在 if 条件语句中。大量的 warp 发散会影响性能，但我没有检查未使用结果的指令是否计入 IPC。

您可能需要研究读取纹理缓存。纹理针对 2D 空间局部性进行了优化，并更好地支持半随机 2D 访问。它可能有助于您的[i][j] 类型访问。

在当前的解决方案中，确保 Y 位置 ([j]) 在具有相邻线程 ID 的两个线程之间变化最小(以尽可能保持内存访问顺序)。

可能是编译器已经为你优化了这个，但是你重新计算了很多次idx*ydim+idy。尝试计算一次并重新使用结果。不过，如果您的算法受计算限制，那将有更大的改进潜力。