超並列プロセサ - GeForceアーキテクチャとCUDAプログラミング

このループの中で、_shared_ float As[BLOCK_SIZE][BLOCK_SIZE];でシェアードメモリ上に部分行列Asの領域を確保する。また、Bsに対しても同様にシェアードメモリ上に領域を確保する。そして、As[ty][tx] = A[a + wA * ty + tx];とBs[ty][tx] = B[b + wB * ty + tx];で、各スレッドは自分が担当する要素のデータをデバイスメモリからシェアードメモリに転送する。さらに両方の転送が終わったことを__syncthreads();で確認してから、次のループでBLOCK_SIZEの長さのベクトルの内積を計算してCsubに足しこんで行く。このように処理を進めて、全てのブロックを処理し終わると、結果をCsubからデバイスメモリのCに格納する。

    // Loop over all the sub-matrices of A and B
    // required to compute the block sub-matrix
    for (int a = aBegin, b = bBegin;
             a <= aEnd;
             a += aStep, b += bStep) {

        // Declaration of the shared memory array As used to
        // store the sub-matrix of A
        __shared__ float As[BLOCK_SIZE][BLOCK_SIZE];

        // Declaration of the shared memory array Bs used to
        // store the sub-matrix of B
        __shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE];

        // Load the matrices from device memory
        // to shared memory; each thread loads
        // one element of each matrix
        As[ty][tx] = A[a + wA * ty + tx];
        Bs[ty][tx] = B[b + wB * ty + tx];

        // Synchronize to make sure the matrices are loaded
        __syncthreads();

        // Multiply the two matrices together;
        // each thread computes one element
        // of the block sub-matrix
        for (int k = 0; k < BLOCK_SIZE; ++k)
            Csub += AS(ty, k) * BS(k, tx);

        // Synchronize to make sure that the preceding
        // computation is done before loading two new
        // sub-matrices of A and B in the next iteration
        __syncthreads();
    }

    // Write the block sub-matrix to device memory;
    // each thread writes one element
    int c = wB * BLOCK_SIZE * by + BLOCK_SIZE * bx;
    C[c + wB * ty + tx] = Csub;
}

このように、結果となる行列Cの各要素を個別のスレッドで計算しているが、BLOCK_SIZE=16の場合、各スレッドは16回の積和演算を行っており、ブロック全体では4096回の積和演算を行っている。この演算に必要なデバイスメモリからシェアードメモリへのデータ転送は、それぞれ256要素のAsとBsの転送だけであり、演算あたりのデータ転送量が小さい効率の良い計算方法となっている。