3

我想用 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 线程块将 (BS_x+1)*(BS_y+1) 全局内存位置读取到共享内存。我相信这个选择比前一个更好，由于使用了共享内存，虽然并不是所有的访问都真正合并，见

分析我的 CUDA 内核的内存访问合并

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

4

2 回答 2

1

感谢您提供分析信息。

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

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

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

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

于 2012-12-13T18:14:10.937 回答
0

我相信在这种情况下，您的第一个并行解决方案会更好，因为每个线程只读取每个全局内存数组元素一次。因此，将这些数组存储在共享内存中不会给您带来预期的改进。

由于在共享内存中存储数据期间可以更好地合并访问全局内存，它可以加快您的程序，但是恕我直言，如果您使用的是 Compute Capability 2.x，并且共享内存的使用也可以降级，那么这可以通过缓存全局内存访问来平衡银行冲突。

于 2012-12-11T17:31:30.977 回答