4.7 优化应用举例
4.7.3 矩阵转置的优化
case 512:
FUNC(reduce4)<T, 512><<< dimGrid, dimBlock, smemSize >>>(d_idata, d_odata); break;
略//256、128、…、2 case 1:
FUNC(reduce4)<T, 1><<< dimGrid, dimBlock, smemSize >>>(d_idata, d_odata); break;
} break;
略//case 5、6、default 与 case4 类同 }
extern "C"
void FUNC(reduceInt)(int size, int threads, int blocks,
int whichKernel, int *d_idata, int *d_odata) {
FUNC(reduce)<int>(size, threads, blocks, whichKernel, d_idata, d_odata);
}
//同理，有 FUNC(reduceFloat)、FUNC(reduceDouble)等
#endif
4.7.3 矩阵转置的优化
Transpose 是一个矩阵转置程序。本节将通过三个版本的 CUDA 矩阵转置程序介绍如何在 访问 global memory 时避免非合并访问和分区冲突问题（partition conflict，见 3.3.3.3 节），以 及如何利用 shared memory 进行线程间通信，实现效率更高的算法。
首先，需要理解如何对矩阵进行棋盘划分实现转置。
如图 4-19 所示的矩阵按照棋盘划分为多个 3*3 的子块。(0,2)块经转置后位于输出矩阵的 (2,0)位置，并且对块内元素也要进行相应转置。
图 4-19 数组的分块转置 下面给出矩阵转置的主机端代码的主要部分。
//设置执行参数
dim3 grid(size_x / BLOCK_DIM, size_y / BLOCK_DIM, 1);
dim3 threads(BLOCK_DIM, BLOCK_DIM, 1);
//调用 transpose_naïve 或者 transpose 内核函数
//transpose_naive<<< grid, threads >>>(d_odata, d_idata, size_x, size_y);
transpose<<< grid, threads >>>(d_odata, d_idata, size_x, size_y);
4.7.3.1 transpose_naive
//---transpose_naive_kernel.cu---
#define BLOCK_DIM 16
// naïve transpose 在向显存写入结果时不满足合并访问条件，严重影响了性能 __global__ void transpose_naive(float *odata, float* idata, int width, int height) { unsigned int xIndex = blockDim.x * blockIdx.x + threadIdx.x;
unsigned int yIndex = blockDim.y * blockIdx.y + threadIdx.y;
if (xIndex < width && yIndex < height) {
unsigned int index_in = xIndex + width * yIndex;
unsigned int index_out = yIndex + height * xIndex;
odata[index_out] = idata[index_in];
} }
对于 256*4096 的数组，其 grid 如图 4-20 所示，图中每个小矩形代表一个 16*16 的 block，
每个 block，负责处理输入数组对应位置上的一块数据。程序中 xIndex 和 yIndex 是线程在整 个 grid 中的位置，而 index_in 和 index_out 分别是线程从输入矩阵读和向输出矩阵写的数据的 地址。
图 4-20 本例 grid 表示
每个 block 中的一个 half warp 按行连续读入数据，在读入数据时尚满足合并访问条件；
但在向输出数组写入结果时，每个 half warp 写入的数据之间的间隔很大，不满足合并访问 条件。此时，一个 half-warp 的访存请求就会产生 16 次单独的传输，这将大大降低可用的显 存带宽。
4.7.3.2 transpose
注意到每个 16*16 分块中，在输入中处于同一列的数据处于输出矩阵同一行中，有可能 可以用合并访问方式写入显存中。通过引入 shared memory 实现线程间通信，就可以让一个 half-warp 中的线程按照合并访问方式输入一行数据后，再以合并访问方式输出数据。
本例中，首先用合并访问方式将数据从显存读入共享存储器中，经过同步后，每个线程 与和它按对角线对称的线程交换操作的数据，再按照合并访问方式将结果写到显存中。这种划 分方式很好地说明了 CUDA 的两层并行：在同一个 block 中实现需要进行数据交换和通信的 细粒度并行，而在各个 block 间实现不需要进行数据交换的粗粒度并行。
/* 这个 kernel 对全局内存器的所有读写均是合并读写，也避免了共享存储器中的 bank conflict 问题。这比 此前的 naïve 版本要快 11 倍。注意，共享存储器大小是(BLOCK_DIM+1)*BLOCK_DIM。注意每一行都加 了 1，这样才能保证 half-warp 按列访问数组的时候不发生 bank conflict。*/
__global__ void transpose(float *odata, float *idata, int width, int height) {
__shared__ float block[BLOCK_DIM][BLOCK_DIM+1]; //静态分配 sharedMem
//把矩阵块读入共享存储器
unsigned int xIndex = blockIdx.x * BLOCK_DIM + threadIdx.x;
unsigned int yIndex = blockIdx.y * BLOCK_DIM + threadIdx.y;
if((xIndex < width) && (yIndex < height)) { //保证内存访问不会超过边界 unsigned int index_in = yIndex * width + xIndex;
block[threadIdx.y][threadIdx.x] = idata[index_in];
}
__syncthreads();
//将转置后的矩阵写回全局存储器
xIndex = blockIdx.y * BLOCK_DIM + threadIdx.x;
yIndex = blockIdx.x * BLOCK_DIM + threadIdx.y;
if((xIndex < height) && (yIndex < width)) {
unsigned int index_out = yIndex * height + xIndex;
odata[index_out] = block[threadIdx.x][threadIdx.y];
} }
注意到 share memory 中的数组 block 大小被设成了 16×17，而不是 16×16。这样每行中处 于同一列的数据就会被存储在不同的 shared memory bank 中，避免了 bank conflict。
transpose 和 transpose_naive 的效率可以相差一个数量级以上。从上面的例子可以看出，使 用 shared memory 可以大大提高某些场合下程序的运行效率。
4.7.3.3 transposeDiagonal
在前两个版本的转置实现中，都存在分区冲突问题（见 3.3.3.3 节）。例如，在 GTX 280GPU 中存在 8 个分区，每个分区中存储连续的 256Byte 数据。如果输入数组的每行中有 2048 个 float 元素，那么 0～7 号元素处于 partition 0 中，而 8～15 号元素处于 partition 1 中，并依次类推，
直到 partition 7 中的第 504～511 个元素。而第 512～519 号元素又被存放在 partition 0 中。
前两个版本中，每个 block 从输入数组中读取一个 32×32 的子块，如图 4-21 所示，block 0 和 block1 读入的数据位于 partition 0 中，block2 和 block 3 读入的数据位于 partition 1 中，并 依次类推。总的来说，各个 block 对各个 partition 的读访问请求的分布基本均衡。
图 4-21 各 block 读取首地址的 partition 分布
但是在写入结果时，刚才的两个例子就发生了严重的分区冲突问题。如图 4-22 所示，block 0～block 63、block 64～block 127 要写入数据的地址都处于第 0 个 partition 中。虽然各个 block 的执行顺序有一定的随机性，但 ID 相近的 block 总的来说还是会集中在同一段时间内执行。
那么这一段时间就只有一个 partition 对应的存储器控制器会处理来自所有 SM 的访存请求，使 得可用显存带宽降低到理论值的 1/8。
图 4-22 各 block 写回首地址的 partition 分布
那么应该如何设计，实现数据访问在多个 partition 间的负载均衡呢？这里给出了转置的第 三个版本。
__global__ void transposeDiagonal(float *odata, float *idata, int width,int height, int nreps) {
__shared__ float tile[TILE_DIM][TILE_DIM+1];
int blockIdx_y= blockIdx.x;
int blockIdx_x= (blockIdx.x+blockIdx.y)%gridDim.x;
int xIndex = blockIdx_x* TILE_DIM + threadIdx.x;
int yIndex = blockIdx_y* TILE_DIM + threadIdx.y;
int index_in = xIndex + (yIndex)*width;
xIndex = blockIdx_y* TILE_DIM + threadIdx.x;
yIndex = blockIdx_x* TILE_DIM + threadIdx.y;
int index_out = xIndex + (yIndex)*height;
for (int r=0; r < nreps; r++) {
for (int i=0; i<TILE_DIM; i+=BLOCK_ROWS) {
tile[threadIdx.y+i][threadIdx.x] = idata[index_in+i*width];
}
__syncthreads();
for (int i=0; i<TILE_DIM; i+=BLOCK_ROWS) {
odata[index_out+i*height] = tile[threadIdx.x][threadIdx.y+i];
} } }
在这个版本的代码中，重点是
int blockIdx_x= (blockIdx.x+blockIdx.y)%gridDim.x;
其中(blockIdx.x+blockIdx.y) 实现输入数据块的增量，取模控制着循环处理一行。图 4-23 和图 4-24 分别是各 block 读和写子块的分区分布。主对角线上的白色矩形是 block(0,0)～
block(63,0) 负 责 的 子 块 ； 淡 灰 色 矩 形 是 block(0,1) ～ block(63,1) 负 责 的 子 块 ； 灰 色 矩 形 是 block(0,2)～block(63,2)负责的子块；其他依此类推。
图 4-23 各 block 负责读入的子块的分区分布
图 4-24 各 block 负责写回的子块的分区分布
从图 4-23 和图 4-24 中可以看出，在这一个版本中，相邻的 block 处理的是一系列子块位 于与矩阵对角线平行的一条线上，从而使得从显存读写数据时在各个分区上的负载均匀分布，
避免了 partition conflict 问题。