CUDA技術を利用したGPUコンピューティングの実際（後編） ―― FFTを利用した光波の伝播（フレネル回折）をGPUで高速計算
ホスト側のプログラム（リスト1）の変更は1個所のみです．先ほどのKernel1の代わりにKernel2を呼び出すように変更します（リスト1の（5）でカーネル関数名を変更）．
Kernel2のプログラムをリスト4に示します．一見するとKernel1よりプログラムが複雑になり，ループも2重ループになっているため計算速度が遅くなりそうな印象を受けます．
__global__ void Kernel2(float *A, float *B, float *C)
{
//GPUでの行列乗算（シェアード・メモリを使用）
int bx = blockIdx.x;
int by = blockIdx.y;
int tx = threadIdx.x;
int ty = threadIdx.y;
float tmp = 0;
__shared__ float As[BLOCK][BLOCK]; // （1）
__shared__ float Bs[BLOCK][BLOCK];
for (int a = 0, b = 0 ; a < WIDTH; a += BLOCK, b += BLOCK) { // （2）
int a_adr = WIDTH * BLOCK * by + a; // （3）
int b_adr = BLOCK * bx + WIDTH * b;
As[ty][tx] = A[a_adr + WIDTH*ty + tx]; // （4）
Bs[ty][tx] = B[b_adr + WIDTH*ty + tx];
__syncthreads(); // （5）
for (int k = 0; k < BLOCK; k++) { // （6）
tmp += As[ty][k] * Bs[k][tx];
}
__syncthreads(); // （7）
}
int adr = WIDTH * BLOCK * by + BLOCK * bx; // （8）
C[adr + WIDTH * ty + tx] = tmp;
}
リスト4 シェアード・メモリを使用するカーネル関数
シェアード・メモリを使用するためには，変数名の前に __shared__ を付ける必要があります（リスト4の（1））．ここでは，BLOCK×BLOCKサイズの2次元配列AsとBsがシェアード・メモリへと割り当てられます．BLOCKは16としてマクロ定義したので，マルチプロセッサ当たりのシェアード・メモリ使用量は2Kバイト（16×16×4バイト×2配列）となります．
Kernel1は，各スレッドに対応する行ベクトルと列ベクトルを読み出してきて，行列の積を計算しました．一方，ここで紹介するKernel2は，行列Aの部分行列（行幅＝WIDTH，列幅＝BLOCK）と行列Bの部分行列（行幅＝BLOCK，列幅＝WIDTH）を読み出し，シェアード・メモリAsとBsへ格納します．ただしAsとBsのサイズはBLOCK×BLOCKなので，一度には格納できません．そこで，これらの部分行列からBLOCK×BLOCKサイズを切り出し，AsとBsへコピーするためのループがリスト4の（2）となっています．（3），（4）で実際に，行列A，Bの部分行列のアドレス（a_adr，b_adr）を算出し，シェアード・メモリへコピーします．
アドレス計算およびコピーの作業は，ブロック内の各スレッドが担当します．各スレッドはばらばらに並列処理されていきます．シェアード・メモリAs，Bsのすべてに正しい値がコピーされる前に次のステップへ進むスレッドを防ぐため，（5）でスレッドの同期を行います．同期にはCUDAのAPIである__syncthreadを使用します．__syncthread には引き数はありません．あるスレッドが先に __syncthread に到着した場合，同じブロック内で動作しているほかのすべてのスレッドがこの __syncthread に到達するまで待機状態になります．全スレッドが __syncthread に到達したあと，スレッドはそれ以降の処理を再開します．
（6）では，ブロック内の各スレッドがAs，Bsから対応する行ベクトルと列ベクトルを読み出し，並列にベクトル積を計算します．（7）で，各スレッドがベクトル積を計算し終えるのを待ち合わせます．
このブロックが担当する行列Aの部分行列と行列Bの部分行列の計算がすべて完了したあと，（8）で行列Cに計算結果を書き込みます．
このカーネルの実行時間は11.3msとなり，Kernel1より高速化できました．
