3

電磁気2D有限差分時間領域（FDTD）コードをCUDA言語で記述したいと思います。磁場の更新のための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も使用しているため、次の2つの理由でこのソリューションに満足できませんでした。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が高くなるため、2番目のアルゴリズムの方が優れています。それでも、CC 2.0では最大IPCは2.0であるため、2番目のソリューションである1.018の平均は、使用可能な計算能力の半分しか使用されていないことを意味します。通常、これはアルゴリズムがメモリにバインドされていることを意味しますが、カーネル内のほとんどすべてのコードが条件内にあるため、あなたの場合はわかりませんif。かなりの量のワープ発散はパフォーマンスに影響しますが、結果が使用されない命令がIPCにカウントされるかどうかは確認していません。

テクスチャキャッシュを読み取ることを検討することをお勧めします。テクスチャは2D空間局所性に最適化されており、セミランダム2Dアクセスをより適切にサポートします。それはあなたの[i][j]タイプアクセスを助けるかもしれません。

[j]現在のソリューションでは、隣接するスレッドIDを持つ2つのスレッド間で最も変化が少ないのはY位置（）であることを確認してください（メモリアクセスを可能な限りシーケンシャルに保つため）。

コンパイラがこれを最適化した可能性がありますが、idx*ydim+idy何度も再計算します。一度計算して、結果を再利用してみてください。ただし、アルゴリズムが計算にバインドされている場合は、改善の可能性が高くなります。

于 2012-12-13T18:14:10.937 に答える
0

この場合、各スレッドが各グローバルメモリ配列要素を1回だけ読み取るため、最初の並列ソリューションの方が優れていると思います。したがって、これらのアレイを共有メモリに保存しても、期待される改善は得られません。

共有メモリにデータを保存する際のグローバルメモリへのアクセスがより適切に統合されるため、プログラムが高速化される可能性がありますが、Compute Capability 2.xを使用している場合は、グローバルメモリアクセスのキャッシュによってバランスが取れており、共有メモリの使用もダウングレードできます。銀行の紛争に。

于 2012-12-11T17:31:30.977 に答える