https://developer.nvidia.com/blog/efficient-matrix-transpose-cuda-cc/
Matrix Transpose
对单精度值矩阵的转置,该矩阵可以在out-of-place的位置进行操作,即输入和输出是内存中的单独数组。为了表示简单,我们只考虑维数是32的整数倍的方阵。它由几个内核以及执行典型任务的主机代码组成,如主机和设备之间的分配和数据传输、启动和计时以及验证其结果、释放主机和设备内存。
除了执行几种不同的矩阵转置之外,我们还运行简单的矩阵复制内核,因为复制性能指示我们希望矩阵转置实现的性能。对于矩阵复制和转置,相关的性能指标是有效带宽,以GB / s表示,方法是将矩阵的单位为GB的大小除以两倍(一次用于加载矩阵,一次用于存储)除以执行秒数。
本研究中的所有内核均启动32×8线程的块(代码中TILE_DIM = 32,BLOCK_ROWS = 8),并且每个线程块转置(或复制)大小为32×32的图块。使用线程数少于平铺中的元素的线程块有利于矩阵转置,因为每个线程转置四个矩阵元素,因此,大部分索引计算成本将在这些元素上摊销。
本例中的内核使用笛卡尔(x,y)映射而不是行/列映射将线程映射到矩阵元素,以简化CUDA C自动变量组件的含义 :threadIdx.x是水平的和threadIdx.y是垂直的。这种映射取决于程序员。要记住的重要一点是,为了确保内存合并,我们希望将变化最快的组件映射到内存中的连续元素。在Fortran中,连续地址对应于多维数组的第一个索引,并且threadIdx.x和blockIdx.x分别在块和网格内变化最快。
Simple Matrix Copy
Naive Matrix Transpose
在transposeNaive中,从idata读取的内容与复制内核中的内容相同,但是对于我们的1024×1024测试矩阵,对odata的写入在连续线程之间具有1024个元素或4096字节的跨度。这使我们很好地陷入了strided的内存访问图的渐近线,并且我们希望此内核的性能会因此受到影响。 copy和transposeNaive内核的结果证明了这一点。
TransposeNaive内核仅达到复制内核的有效带宽的一小部分。由于此内核除了复制外几乎没有其他作用,因此我们希望更接近复制吞吐量。让我们看看我们如何做到这一点。
Coalesced Transpose Via Shared Memory
较差的转置性能的补救方法是使用共享内存,以避免通过全局内存大步前进。下图描述了转置中如何使用共享内存。
#include <stdio.h>
#include <assert.h>
// Convenience function for checking CUDA runtime API results
// can be wrapped around any runtime API call. No-op in release builds.
inline
cudaError_t checkCuda(cudaError_t result)
{
#if defined(DEBUG) || defined(_DEBUG)
if (result != cudaSuccess) {
fprintf(stderr, "CUDA Runtime Error: %s\n", cudaGetErrorString(result));
assert(result == cudaSuccess);
}
#endif
return result;
}
const int TILE_DIM = 32;
const int BLOCK_ROWS = 8;
const int NUM_REPS = 100;
// Check errors and print GB/s
void postprocess(const float* ref, const float* res, int n, float ms)
{
bool passed = true;
for (int i = 0; i < n; i++)
if (res[i] != ref[i]) {
printf("%d %f %f\n", i, res[i], ref[i]);
printf("%25s\n", "*** FAILED ***");
passed = false;
break;
}
if (passed)
printf("%20.2f\n", 2 * n * sizeof(float) * 1e-6 * NUM_REPS / ms);
}
// simple copy kernel
// Used as reference case representing best effective bandwidth.
__global__ void copy(float* odata, const float* idata)
{//数据拷贝
int x = blockIdx.x * TILE_DIM + threadIdx.x;
int y = blockIdx.y * TILE_DIM + threadIdx.y;//还请注意,必须在矩阵索引y的计算中使用TILE_DIM,而不是BLOCK_ROWS或blockDim.y。
int width = gridDim.x * TILE_DIM;
for (int j = 0; j < TILE_DIM; j += BLOCK_ROWS)//每个线程在此例程结束时在循环中复制矩阵的四个元素,因为块中的线程数比tile中的元素数小四倍(TILE_DIM / BLOCK_ROWS)。
odata[(y + j) * width + x] = idata[(y + j) * width + x];// 循环在第二个维度而不是第一个维度上进行迭代,以便连续线程加载和存储连续数据,并且对idata的所有读取和对odata的写入均被合并。
}
// copy kernel using shared memory
// Also used as reference case, demonstrating effect of using shared memory.
__global__ void copySharedMem(float* odata, const float* idata)
{
__shared__ float tile[TILE_DIM * TILE_DIM];
int x = blockIdx.x * TILE_DIM + threadIdx.x;
int y = blockIdx.y * TILE_DIM + threadIdx.y;
int width = gridDim.x * TILE_DIM;
for (int j = 0; j < TILE_DIM; j += BLOCK_ROWS)
tile[(threadIdx.y + j) * TILE_DIM + threadIdx.x] = idata[(y + j) * width + x];
__syncthreads();
for (int j = 0; j < TILE_DIM; j += BLOCK_ROWS)
odata[(y + j) * width + x] = tile[(threadIdx.y + j) * TILE_DIM + threadIdx.x];
}
// naive transpose
// Simplest transpose; doesn't use shared memory.
// Global memory reads are coalesced but writes are not.
__global__ void transposeNaive(float* odata, const float* idata)
{
int x = blockIdx.x * TILE_DIM + threadIdx.x;
int y = blockIdx.y * TILE_DIM + threadIdx.y;
int width = gridDim.x * TILE_DIM;
for (int j = 0; j < TILE_DIM; j += BLOCK_ROWS)
odata[x * width + (y + j)] = idata[(y + j) * width + x];//我们的第一个转置内核看起来与复制内核非常相似。唯一的区别是交换了odata的索引。
}
// coalesced transpose
// Uses shared memory to achieve coalesing in both reads and writes
// Tile width == #banks causes shared memory bank conflicts.
__global__ void transposeCoalesced(float* odata, const float* idata)
{
// 以下内核执行此“tiled”转置。transposeCoalesced结果是对transposeNaive情况的改进,但与copy内核的性能仍然相去甚远。//我们可能会猜测,造成性能差距的原因是与使用共享内存和所需的同步阻塞__syncthreads()相关的开销。
__shared__ float tile[TILE_DIM][TILE_DIM];//我们可以使用以下使用共享内存的复制内核轻松地对此进行测试。
//对于32×32元素的共享内存tile,数据列中的所有元素都映射到同一共享内存库,从而导致内存库冲突的最坏情况:读取数据列将导致32路存储库冲突。
int x = blockIdx.x * TILE_DIM + threadIdx.x;
int y = blockIdx.y * TILE_DIM + threadIdx.y;
int width = gridDim.x * TILE_DIM;
for (int j = 0; j < TILE_DIM; j += BLOCK_ROWS)
tile[threadIdx.y + j][threadIdx.x] = idata[(y + j) * width + x];//在第一个do循环中,线程warp将idata中的连续数据读取到共享内存tile的行中。
__syncthreads();//因为线程向odata写入的数据与从idata读取的数据不同,所以我们必须使用逐块阻塞同步__syncthreads()。
//请注意,在这种情况下,从技术上讲不需要synchthreads()调用,因为对元素的操作是由同一线程执行的,但是我们在此处包括了它,以模仿转置行为。
x = blockIdx.y * TILE_DIM + threadIdx.x; // transpose block offset
y = blockIdx.x * TILE_DIM + threadIdx.y;
for (int j = 0; j < TILE_DIM; j += BLOCK_ROWS)
odata[(y + j) * width + x] = tile[threadIdx.x][threadIdx.y + j];
}
// No bank-conflict transpose
// Same as transposeCoalesced except the first tile dimension is padded
// to avoid shared memory bank conflicts.
__global__ void transposeNoBankConflicts(float* odata, const float* idata)
{
__shared__ float tile[TILE_DIM][TILE_DIM + 1];
int x = blockIdx.x * TILE_DIM + threadIdx.x;
int y = blockIdx.y * TILE_DIM + threadIdx.y;
int width = gridDim.x * TILE_DIM;
for (int j = 0; j < TILE_DIM; j += BLOCK_ROWS)
tile[threadIdx.y + j][threadIdx.x] = idata[(y + j) * width + x];
__syncthreads();
x = blockIdx.y * TILE_DIM + threadIdx.x; // transpose block offset
y = blockIdx.x * TILE_DIM + threadIdx.y;
for (int j = 0; j < TILE_DIM; j += BLOCK_ROWS)
odata[(y + j) * width + x] = tile[threadIdx.x][threadIdx.y + j];
}
int main(int argc, char** argv)
{
const int nx = 1024;
const int ny = 1024;
const int mem_size = nx * ny * sizeof(float);
dim3 dimGrid(nx / TILE_DIM, ny / TILE_DIM, 1);
dim3 dimBlock(TILE_DIM, BLOCK_ROWS, 1);
int devId = 0;
if (argc > 1) devId = atoi(argv[1]);
cudaDeviceProp prop;
checkCuda(cudaGetDeviceProperties(&prop, devId));
printf("\nDevice : %s\n", prop.name);
printf("Matrix size: %d %d, Block size: %d %d, Tile size: %d %d\n",
nx, ny, TILE_DIM, BLOCK_ROWS, TILE_DIM, TILE_DIM);
printf("dimGrid: %d %d %d. dimBlock: %d %d %d\n",
dimGrid.x, dimGrid.y, dimGrid.z, dimBlock.x, dimBlock.y, dimBlock.z);
checkCuda(cudaSetDevice(devId));
float* h_idata = (float*)malloc(mem_size);
float* h_cdata = (float*)malloc(mem_size);
float* h_tdata = (float*)malloc(mem_size);
float* gold = (float*)malloc(mem_size);
float* d_idata, * d_cdata, * d_tdata;
checkCuda(cudaMalloc(&d_idata, mem_size));
checkCuda(cudaMalloc(&d_cdata, mem_size));
checkCuda(cudaMalloc(&d_tdata, mem_size));
// check parameters and calculate execution configuration
if (nx % TILE_DIM || ny % TILE_DIM) {
printf("nx and ny must be a multiple of TILE_DIM\n");
goto error_exit;
}
if (TILE_DIM % BLOCK_ROWS) {
printf("TILE_DIM must be a multiple of BLOCK_ROWS\n");
goto error_exit;
}
// host
for (int j = 0; j < ny; j++)
for (int i = 0; i < nx; i++)
h_idata[j * nx + i] = j * nx + i;
// correct result for error checking
for (int j = 0; j < ny; j++)
for (int i = 0; i < nx; i++)
gold[j * nx + i] = h_idata[i * nx + j];
// device
checkCuda(cudaMemcpy(d_idata, h_idata, mem_size, cudaMemcpyHostToDevice));
// events for timing
cudaEvent_t startEvent, stopEvent;
checkCuda(cudaEventCreate(&startEvent));
checkCuda(cudaEventCreate(&stopEvent));
float ms;
// ------------
// time kernels
// ------------
printf("%25s%25s\n", "Routine", "Bandwidth (GB/s)");
// ----
// copy
// ----
printf("%25s", "copy");
checkCuda(cudaMemset(d_cdata, 0, mem_size));
// warm up
copy << <dimGrid, dimBlock >> > (d_cdata, d_idata);
checkCuda(cudaEventRecord(startEvent, 0));
for (int i = 0; i < NUM_REPS; i++)
copy << <dimGrid, dimBlock >> > (d_cdata, d_idata);
checkCuda(cudaEventRecord(stopEvent, 0));
checkCuda(cudaEventSynchronize(stopEvent));
checkCuda(cudaEventElapsedTime(&ms, startEvent, stopEvent));
checkCuda(cudaMemcpy(h_cdata, d_cdata, mem_size, cudaMemcpyDeviceToHost));
postprocess(h_idata, h_cdata, nx * ny, ms);
// -------------
// copySharedMem
// -------------
printf("%25s", "shared memory copy");
checkCuda(cudaMemset(d_cdata, 0, mem_size));
// warm up
copySharedMem << <dimGrid, dimBlock >> > (d_cdata, d_idata);
checkCuda(cudaEventRecord(startEvent, 0));
for (int i = 0; i < NUM_REPS; i++)
copySharedMem << <dimGrid, dimBlock >> > (d_cdata, d_idata);
checkCuda(cudaEventRecord(stopEvent, 0));
checkCuda(cudaEventSynchronize(stopEvent));
checkCuda(cudaEventElapsedTime(&ms, startEvent, stopEvent));
checkCuda(cudaMemcpy(h_cdata, d_cdata, mem_size, cudaMemcpyDeviceToHost));
postprocess(h_idata, h_cdata, nx * ny, ms);
// --------------
// transposeNaive
// --------------
printf("%25s", "naive transpose");
checkCuda(cudaMemset(d_tdata, 0, mem_size));
// warmup
transposeNaive << <dimGrid, dimBlock >> > (d_tdata, d_idata);
checkCuda(cudaEventRecord(startEvent, 0));
for (int i = 0; i < NUM_REPS; i++)
transposeNaive << <dimGrid, dimBlock >> > (d_tdata, d_idata);
checkCuda(cudaEventRecord(stopEvent, 0));
checkCuda(cudaEventSynchronize(stopEvent));
checkCuda(cudaEventElapsedTime(&ms, startEvent, stopEvent));
checkCuda(cudaMemcpy(h_tdata, d_tdata, mem_size, cudaMemcpyDeviceToHost));
postprocess(gold, h_tdata, nx * ny, ms);
// ------------------
// transposeCoalesced
// ------------------
printf("%25s", "coalesced transpose");
checkCuda(cudaMemset(d_tdata, 0, mem_size));
// warmup
transposeCoalesced << <dimGrid, dimBlock >> > (d_tdata, d_idata);
checkCuda(cudaEventRecord(startEvent, 0));
for (int i = 0; i < NUM_REPS; i++)
transposeCoalesced << <dimGrid, dimBlock >> > (d_tdata, d_idata);
checkCuda(cudaEventRecord(stopEvent, 0));
checkCuda(cudaEventSynchronize(stopEvent));
checkCuda(cudaEventElapsedTime(&ms, startEvent, stopEvent));
checkCuda(cudaMemcpy(h_tdata, d_tdata, mem_size, cudaMemcpyDeviceToHost));
postprocess(gold, h_tdata, nx * ny, ms);
// ------------------------
// transposeNoBankConflicts
// ------------------------
printf("%25s", "conflict-free transpose");
checkCuda(cudaMemset(d_tdata, 0, mem_size));
// warmup
transposeNoBankConflicts << <dimGrid, dimBlock >> > (d_tdata, d_idata);
checkCuda(cudaEventRecord(startEvent, 0));
for (int i = 0; i < NUM_REPS; i++)
transposeNoBankConflicts << <dimGrid, dimBlock >> > (d_tdata, d_idata);
checkCuda(cudaEventRecord(stopEvent, 0));
checkCuda(cudaEventSynchronize(stopEvent));
checkCuda(cudaEventElapsedTime(&ms, startEvent, stopEvent));
checkCuda(cudaMemcpy(h_tdata, d_tdata, mem_size, cudaMemcpyDeviceToHost));
postprocess(gold, h_tdata, nx * ny, ms);
error_exit:
// cleanup
checkCuda(cudaEventDestroy(startEvent));
checkCuda(cudaEventDestroy(stopEvent));
checkCuda(cudaFree(d_tdata));
checkCuda(cudaFree(d_cdata));
checkCuda(cudaFree(d_idata));
free(h_idata);
free(h_tdata);
free(h_cdata);
free(gold);
}
在这篇文章中,我们介绍了三个代表内核转置的各种优化的内核。内核展示了如何使用共享内存来合并全局内存访问,以及如何填充阵列以避免共享内存库冲突。从我们的内核的相对收益来看,coalescing global memory访问是迄今为止获得良好性能的最关键方面,这在许多应用程序中都是如此。由于全局内存合并非常重要,因此当我们在3D网格上查看有限差分计算时,我们将在下一篇文章中再次进行讨论。