シェアードメモリとバンク

GPUを使った行列転置でほぼ頻出話題のバンクコンフリクト(Bank Conflict)とパディング(Padding)。前提知識としてCUDA等のGPUプログラミングでシェアードメモリ(共有メモリ)を明示的に扱う際に、Bankを意識しないと遅くなっちゃいますよという話があります。
bank

シェアードメモリは4byte刻みで32個のBankに分かれています。
例えばthread 0がSMEM[0]、thread 1がSMEM[32]にアクセスするようなプログラムだと(SMEMは4byte型とします)、どちらも同じBank0にアクセスすることになるのでリプレイ(replay)が発生して遅くなります。

行列転置バンクコンフリクトありver

つまり普通にシェアードメモリを使って行列転置をしようとすると
jurai

こんな感じにShared Memoryから値を読み込む際にバンクコンフリクトが起こってしまっています。
ちなみにここではわかりやすくBankは0～3までとしています。

行列転置パディングありver

そこでパディングをいれることで参照Bankがずれるのでこれを解消できます。そこそこ有名なテクニックです。
Kepler GPUアーキテクチャとプログラム最適化
プログラムを高速化する話Ⅱ 〜GPGPU編〜

sample0

概念図
juraip

これでShared Memoryに書き込むときも読み込むときも、Bank競合が起きてないことが分かります。
でも貴重なShared Memoryに未使用領域があるのはちょっと気持ち悪くないですか。

行列転置パディングなしver

上の概念図がすでに大きなヒントでした。
sinki
これでも読み書きともバンクコンフリクトは起こってないですね。
メモリアクセスの仕方としては
code0

こんな感じにx方向のズレ幅をyごとにずらしてやればいいです。

一応コードの全貌を載せておきます。

#include<stdio.h>
#include<cuda.h>
#include<cuda_runtime.h>

// --matrix size --
const int nimax=8192,njmax=8192;
// -- the number of threads per block --
const int BS=32;
__device__ float idata[nimax][njmax],odata[nimax][njmax];

__global__ void transposeNaive(){
  int ni,nj;
  ni=blockIdx.y*blockDim.y+threadIdx.y;
  nj=blockIdx.x*blockDim.x+threadIdx.x;
  odata[nj][ni]=idata[ni][nj];
  return;
}

__global__ void transposeShared(){
  int ni,nj;
  __shared__ float tile[BS][BS];
  ni=blockIdx.y*blockDim.y+threadIdx.y;
  nj=blockIdx.x*blockDim.x+threadIdx.x;
  tile[threadIdx.y][threadIdx.x]=idata[ni][nj];
  __syncthreads();
  ni=blockIdx.x*blockDim.x+threadIdx.y;
  nj=blockIdx.y*blockDim.y+threadIdx.x;
  odata[ni][nj]=tile[threadIdx.x][threadIdx.y];
  return;
}

__global__ void transposeSharedNobankconf(){
  int ni,nj;
  __shared__ float tile[BS][BS+1];
  ni=blockIdx.y*blockDim.y+threadIdx.y;
  nj=blockIdx.x*blockDim.x+threadIdx.x;
  tile[threadIdx.y][threadIdx.x]=idata[ni][nj];
  __syncthreads();
  ni=blockIdx.x*blockDim.x+threadIdx.y;
  nj=blockIdx.y*blockDim.y+threadIdx.x;
  odata[ni][nj]=tile[threadIdx.x][threadIdx.y];
  return;
}


__global__ void transposeSharedNobankconf_nopading(){
  int ni,nj;
  __shared__ float tile[BS][BS];
  ni=blockIdx.y*blockDim.y+threadIdx.y;
  nj=blockIdx.x*blockDim.x+threadIdx.x;
  tile[threadIdx.y][(threadIdx.x+threadIdx.y)%BS]=idata[ni][nj];
  __syncthreads();
  ni=blockIdx.x*blockDim.x+threadIdx.y;
  nj=blockIdx.y*blockDim.y+threadIdx.x;
  odata[ni][nj]=tile[threadIdx.x][(threadIdx.x+threadIdx.y)%BS];
  return;
}


void transpose(float hidata[][njmax],float hodata[][njmax]){
  int ni,nj;
  for(ni=0;ni<nimax;ni++){
    for(nj=0;nj<njmax;nj++){
      hodata[nj][ni]=hidata[ni][nj];
    }
  }
  return;
}

void check_mat(float C[][njmax],float CC[][njmax]){
  int ni,nj,flg;
  flg=0;
  for(ni=0;ni<nimax;ni++){
    for(nj=0;nj<njmax;nj++){
      if(fabs(C[ni][nj]-CC[ni][nj])>1.0e-8){
        flg=1;
        printf("%d %d %lf %lf\n",ni,nj,C[ni][nj],CC[ni][nj]);
      }
    }
  }
  if(flg==1){
    printf("Calculation error.\n");
  }else{
    printf("OK.\n");
  }
  return;
}




float hidata[nimax][njmax],hodata[nimax][njmax],htdata[nimax][njmax];

int main(){
  int ni,nj;
  dim3 grid,blck;
  grid.x=nimax/BS; grid.y=njmax/BS;
  blck.x=BS; blck.y=BS;
  // -- set initial value --
  srand(248309);
  
  
  for(ni=0;ni<nimax;ni++)
  {
    for(nj=0;nj<njmax;nj++)
    {
      hidata[ni][nj]=(float)rand()*1.1;
    }
  }
  
  
  cudaMemcpyToSymbol(idata,hidata,nimax*njmax*sizeof(float));
  transpose(hidata,htdata);
  // -- 1. transpose by simple kernel --
  transposeNaive<<<grid,blck>>>();
  transposeNaive<<<grid,blck>>>();
  cudaMemcpyFromSymbol(hodata,odata,nimax*njmax*sizeof(float));
  check_mat(htdata,hodata);
  
  // -- 2. transpose with shared memory --
  transposeShared<<<grid,blck>>>();
  transposeShared<<<grid,blck>>>();
  cudaMemcpyFromSymbol(hodata,odata,nimax*njmax*sizeof(float));
  check_mat(htdata,hodata);
  
  // -- 3. transpose with shared memory and without bank conflict --
  transposeSharedNobankconf<<<grid,blck>>>();
  transposeSharedNobankconf<<<grid,blck>>>();
  cudaMemcpyFromSymbol(hodata,odata,nimax*njmax*sizeof(float));
  check_mat(htdata,hodata);

  // -- 4. transpose with shared memory and without bank conflict and nopading --
  transposeSharedNobankconf_nopading<<<grid,blck>>>();
  transposeSharedNobankconf_nopading<<<grid,blck>>>();
  cudaMemcpyFromSymbol(hodata,odata,nimax*njmax*sizeof(float));
  check_mat(htdata,hodata);
  

  cudaDeviceReset();
  
  return 0;
}

8192×8192要素の行列の転置をするコードとなっています。

速度比較


normal

1.transposeNaive：シェアードメモリすら使ってない。グローバルメモリへコアレスアクセスできてないので一番遅い
2.transposeShared：シェアードメモリ内でバンクコンフリクト発生
3.transposeNobankconf：パディング入れてるやつ
4.transposeNobankconf_nopading：パディング入れてないやつ


まずグローバルメモリにコアレスアクセスできてない1番、これは論外です。次に2,3のバンクコンフリクトのあるなし、これだけで2倍違うことがわかります。
今回考案した4番は、速度的に3番に劣ってないことがわかります。

で、これは何の役に立つんですか？

さきほど書いたよう、パディング入りは無駄なメモリ領域を確保しています。なのでOccupancyが100%になる前にシェアードメモリ容量が律速になるような問題では以下のようなことが発生します。

1Blockあたり16.25KBのシェアードメモリを使うプログラムの場合
nvvp2

1Blockあたり16.0KBのシェアードメモリを使うプログラムの場合
nvvp

このようにOccupancyが上がることが期待できます。

なお、上記の行列転置のコードのカーネルはすべてOccupancy100%です。(Shared Memoryを使い切る前にOccupancy100%になる)