https://dmacssite.github.io/materials/MatrixTranspose.pdf
https://dmacssite.github.io/materials/MatrixTranspose.pdf
Chapter 1. Introduction
矩阵转置优化CUDA内存管理
本文档讨论了CUDA应用程序性能的各个方面，这些方面与有效使用GPU内存和应用于矩阵转置的数据管理有关。特别地，本文档讨论了以下内存使用问题:
- 合并数据传输到和从全局内存
- 共享内存库冲突
- 分区冲突
主机和设备之间的数据传输，以及常量和纹理存储器。这里没有讨论高效内存使用的其他方面，例如合并和分区冲突都处理全局设备和片上内存之间的数据传输，而共享内存库冲突处理片上共享内存。
读者应该熟悉基本的CUDA编程概念，如内核、线程和块，以及对CUDA线程可访问的不同内存空间的基本理解。CUDA编程指南以及CUDAZone(http://www.nvidia.com/cuda)上的其他资源提供了对CUDA编程的很好的介绍。
接下来给出矩阵转置问题陈述，然后简要讨论性能指标，之后文档的其余部分呈现一系列CUDA矩阵转置内核，逐步解决各种性能瓶颈。
矩阵转置特性
在本文档中，我们优化了浮点数矩阵的转置操作，即输入和输出矩阵分别位于不同的内存位置。为了表示的简单性和简洁性，我们只考虑其维度为32的整数倍的方阵，即切片大小，通过文档。然而，修改代码以适应任意大小的矩阵是很简单的。
代码突出显示和性能度量
所有转置情况的主机代码在附录a中给出。主机代码执行典型的任务:主机和设备之间的数据分配和传输，几个内核的启动和定时，结果验证，以及主机和设备内存的释放。
除了不同的矩阵转置，我们还运行执行矩阵复制的内核。矩阵副本的性能作为我们希望矩阵转置达到的基准。
对于矩阵复制和转置，相关的性能指标是有效带宽，以GB/s为单位计算为矩阵大小的两倍---次用于读取矩阵，一次用于写入矩阵-一除以执行时间。由于计时是在执行NUM REPS次数的循环中执行的，这是在代码顶部定义的，因此有效带宽也由NUM_REPS规范化。
在代码上循环num rep时间以进行测量有两种不同的方式:在内核启动上循环，以及在内核内循环加载和存储。这些测量的主机代码如下:
cudaEventRecord(start, 0); for (int i=0; i < NUM_REPS; i++) { kernel<<<grid, threads>>>(d_odata, d_idata,size_x,size_y,1); } cudaEventRecord(stop, 0); cudaEventSynchronize(stop); float outerTime; cudaEventElapsedTime(&outerTime, start, stop); ... // take measurements for loop inside kernel cudaEventRecord(start, 0); kernel<<<grid,threads>>> (d_odata, d_idata, size_x, size_y, NUM_REPS); cudaEventRecord(stop, 0); cudaEventSynchronize(stop); float innerTime; cudaEventElapsedTime(&innerTime, start, stop);第一次计时是通过主机代码中的for循环完成的，第二次计时是通过将NUM REPS作为参数传递给内核完成的。一个简单的复制内核如下所示:
__global__ void copy(float *odata, float* idata, int width, int height, int nreps)
{ int xIndex = blockIdx.x*TILE_DIM + threadIdx.x; int yIndex = blockIdx.y*TILE_DIM + threadIdx.y; int index = xIndex + width*yIndex; for (int r=0; r < nreps; r++) { for (int i=0; i<TILE_DIM; i+=BLOCK_ROWS) { odata[index+i*width] = idata[index+i*width]; } }
}这两种时间的区别在于内核启动的开销，这在不同的内核之间应该是一致的，以及在每个内核开始时计算矩阵索引所花费的时间。此外，在内核启动上循环也可以作为一种同步机制。当内核从宿主代码中的循环中多次启动时，一个内核启动的所有块必须在下一次启动的任何块开始之前完成执行。因此，每次循环选代都会重置活动块集和内存访问模式。当在内核内执行循环时，活动线程块集在执行定时循环过程中有更多的机会分散。
计时代码的两种方法都提供了有用的度量方法，第一种方法指示通常使用什么作为总体性能度量，第二种方法用于比较内核之间的数据移动时间。
在下一节中，我们将介绍从主机代码调用的不同内核，每个内核都解决不同的性能问题。本研究中的所有内核都启动尺寸为32x8的线程块，其中每个块转置(或复制)尺寸为32x32的块。因此，参数TILE DIM和BLOCKROWS分别设置为32和8。使用线程数少于tile中元素数的线程块对于矩阵转置是有利的，因为每个线程转置几个矩阵元素，在我们的示例中是四个，并且计算索引的大部分成本是在这些元素上平摊的。
2.复制与转置内核
简单的复制
我们考虑的前两种情况是naive转置和简单复制，每种情况都在32x32矩阵切片上使用32x8线程块。前一节给出了复制内核，它显示了所有内核的基本布局。前两个参数odata和data是指向输入和输出矩阵的指针，width和height是矩阵x和y的维度，nreps决定在矩阵之间执行数据移动的循环次数。在这个内核中，计算全局2D矩阵索引xIndex和yIndex，它们依次用于计算index，即每个线程访问矩阵元素所使用的1D索引。i上的循环为index添加了额外的偏移量，以便每个线程复制数组的多个元素，r上的循环用于多次计时数据从输入到输出数组的传输。
__global__ void transposeNaive(float *odata, float* idata, int width, int height, int nreps)
{ int xIndex = blockIdx.x*TILE_DIM + threadIdx.x; int yIndex = blockIdx.y*TILE_DIM + threadIdx.y; int index_in = xIndex + width * yIndex; int index_out = yIndex + height * xIndex; for (int r=0; r < nreps; r++) { for (int i=0; i<TILE_DIM; i+=BLOCK_ROWS) { odata[index_out+i] = idata[index_in+i*width];} }
}Naïve transpose
The naïve transpose:
__global__ void transposeNaive(float *odata, float* idata, int width, int height, int nreps)
{ int xIndex = blockIdx.x*TILE_DIM + threadIdx.x; int yIndex = blockIdx.y*TILE_DIM + threadIdx.y; int index_in = xIndex + width * yIndex; int index_out = yIndex + height * xIndex; for (int r=0; r < nreps; r++) { for (int i=0; i<TILE_DIM; i+=BLOCK_ROWS) { odata[index_out+i] = idata[index_in+i*width];} }
}几乎与上面的复制内核相同，只是index(用于访问复制内核的输入和输出数组中的元素的数组索引)被两个索引index in(相当于复制内核中的index)和index_out取代。每个执行内核的线程将四个元素从输入矩阵的一列转置到它们在输出矩阵的一行中的转置位置。
这两个内核在2048 x2048矩阵上使用GTX280的性能如下表所示:
复制核和原始转置核之间代码的微小差异对性能有深远的影响--几乎有两个数量级的影响。这就引出了我们的第一个优化技术:全局内存合并。
合并转置
由于设备内存比片上内存具有更高的延迟和更低的带宽，因此必须特别注意如何执行全局内存访问，在我们的示例中，从data加载数据并将数据存储在odata中。如果满足某些条件，由半曲线程进行的所有全局内存访问可以合并到一个或两个事务中。这些标准取决于设备的计算能力，这可以通过运行deviceQuerySDK示例来确定。对于1.0和1.1的计算能力，合并需要满足以下条件:
- 线程必须访问32-64位或128位字，导致一个事务(用于32位和64位字)或两个事务(用于128位字)
- 对于32位和64位字，所有16个字必须位于相同的64字节或128字节的对齐段中，对于128位字，数据必须位于两个连续的128字节对齐段中
- 线程需要按顺序访问单词。如果第k个线程要访问一个单词，那么它必须访问第k个单词，尽管并非所有线程都需要参与。
对于计算能力为1.2的设备，对合并的要求比较宽松。当数据位于32、64和128字节对齐的段中时，无论段内线程的访问模式如何，都可以合并到单个事务中。通常，如果一半的线程访问N个内存段，则发出N个内存事务。
简而言之，如果内存访问合并到计算能力为1.0或1.1的设备上，那么它将合并到计算能力为1.2或更高的设备上。如果它不能在具有1.0或1.1计算能力的设备上合并，那么它可能会合并在计算能力为1.2或更高的设备上，要么完全合并，要么可能导致内存事务数量减少。
对于简单复制和naive转置，来自数据的所有负载都合并到具有上述任何计算能力的设备上。对于i循环中的每次迭代，每次半warp读取16个连续的32位单词，或者读取tile的一半行。通过cudaMalloc()分配设备内存，并选择TILEDIM为16的倍数，确保与内存段对齐，因此所有负载都被合并。
当写入odata时，合并行为在简单复制和naive转置内核之间是不同的。对于简单的复制，在illoop的每次迭代期间，halfwarp以合并的方式写入tile的一半行。在naive转置的情况下，对于i循环的每次迭代，halfwarp将一列浮点数的一半写入不同的内存段，从而产生16个独立的内存事务，而不管计算能力如何。
避免非合并全局内存访问的方法是将数据读入共享内存，并让每个半曲访问共享内存中的不连续位置，以便将连续数据写入odata。共享内存中的不连续访问式不像在全局内存中那样有性能损失，但是上面的过程要求内存中的每个元素由不同的线程访问，因此需要调用a_synchthreads()来确保从数据到共享内存的所有读取都在从共享内存到odata的写入开始之前完成。合并转置列如下:
__global__ void transposeCoalesced(float *odata, float *idata, int width, int height, int nreps)
{ __shared__ float tile[TILE_DIM][TILE_DIM]; 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]; } }
}下面给出了合并转置核中半翘曲的数据流的描述。将数据矩阵瓦片的四个半行写入共享内存32x32数组“切片”由黄线段表示。在调用a_syncthreads()以确保对tile的所有写入都完成之后，halfwarp将tile的四半列写入odata矩阵tile的四半行，由绿色线段表示。
通过改进odata对内存的访问模式，写操作被合并，我们看到了性能的提高
虽然合并转置的有效带宽比原始转置显著增加，但合并转置与副本之间仍然存在很大的性能差距。转置所需的额外索引似乎并不是造成性能差距的原因，因为“内核循环”列中的结果也显示了很大的性能差异，其中索引计算是在数据移动的100次迭代中平摊的。造成这种性能差距的一个可能原因是合并转置中所需的同步屏障。使用下面的复制内核可以很容易地评估这一点，它利用共享内存并包含一个
Syncthreads()调用:
__global__ void copySharedMem(float *odata, float *idata, int width, int height, int nreps)
{ __shared__ float tile[TILE_DIM][TILE_DIM]; int xIndex = blockIdx.x*TILE_DIM + threadIdx.x; int yIndex = blockIdx.y*TILE_DIM + threadIdx.y; int index = xIndex + width*yIndex; 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+i*width]; } __syncthreads(); for (int i=0; i<TILE_DIM; i+=BLOCK_ROWS) { odata[index+i*width] = tile[threadIdx.y+i][threadIdx.x]; } }
} 这个内核的成功执行不需要_syncthreads()调用，因为线程不共享数据，并且只用于评估合并转置中同步屏障的成本。结果如下表所示:
共享内存复制结果似乎表明，在比较简单复制和共享内存复制时，“内核循环”列表明，使用带有同步屏障的共享内存对性能几乎没有影响。然而，在比较合并转置和共享内存复制内核时，关于如何访问共享内存有一个需要解决的性能瓶颈:共享内存库冲突。
共享内存库冲突
共享内存被分成16个大小相等的内存模块，称为存储库，这些存储库的组织方式是将连续的32位字分配给连续的存储库。这些银行可以同时被访问，为了获得最大的带宽进出共享内存，半曲的线程应该访问与不同银行相关联的共享内存。此规则的例外情况是，当半warp中的所有线程读取相同的共享内存地址时，这会导致广播，其中该地址的数据在一个事务中发送给半warp的所有线程。
在分析CUDA应用程序时，可以使用warpserialize标志来确定共享内存库冲突是否发生在任何内核中。一般来说，这个标志也反映了原子和常量内存的使用，但是在我们的示例中这两者都不存在。
合并转置使用32x32的浮点数共享内存数组。对于这个大小的数组，列k和k+16中的所有数据都映射到同一个库。因此，当从共享内存中的tile写入部分列到odata中的行时半warp会经历16路银行冲突并序列化请求。避免这种冲突的一个简单方法是将共享内存数组填充一列:
__shared__ float tile[TILE_DIM][TILE_DIM+1];当向共享内存写入半曲时，填充不会影响共享内存库访问模式，这仍然没有冲突，但是现在通过添加单个列，对列中半曲数据的访问也没有冲突。内核的性能，现在合并和内存库冲突无，添加到我们的下表:
虽然填充共享内存数组确实消除了共享内存库冲突，正如用CUDA分析器检查warp serialize标志所证实的那样，但它对性能的影响很小(在这个阶段实现时)。因此，合并和共享内存库的无冲突转置与共享内存内存复制之间仍然存在很大的性能差距。在下一节中，我们将把转置分解为多个组件，以确定导致性能下降的原因。
分解的转置
在最佳优化的转置和上表中的共享内存副本之间存在超过4倍的性能差异。这种情况不仅适用于在内核启动时进行循环的测量，也适用于在内核内进行循环的测量，其中与附加索引计算相关的成本在100次迭代中平摊。
为了进一步研究，我们重新审视转置的数据流，并将其与副本的数据流进行比较，这两者都在下面的图表的顶部表示。复制代码和转置代码本质上有两个不同之处:将数据转置到一个tile内，并将数据写入转置的tile中。我们可以通过实现两个单独执行其中一个组件的内核来隔离这两个组件之间的性能。如下图的下半部分所示，细粒度转置内核将数据转置到一个块内，但将该块写入副本将写入该块的位置。粗粒度转置内核将转置块写入odata矩阵中的转置位置，但不转置块内的数据。
这两个内核的源代码如下:
__global__ void transposeFineGrained(float *odata, float *idata, int width, int height, int nreps)
{ __shared__ float block[TILE_DIM][TILE_DIM+1]; int xIndex = blockIdx.x * TILE_DIM + threadIdx.x;int yIndex = blockIdx.y * TILE_DIM + threadIdx.y;int index = xIndex + (yIndex)*width; for (int r=0; r<nreps; r++) { for (int i=0; i < TILE_DIM; i += BLOCK_ROWS) { block[threadIdx.y+i][threadIdx.x] = idata[index+i*width]; } __syncthreads(); for (int i=0; i < TILE_DIM; i += BLOCK_ROWS) { odata[index+i*height] = block[threadIdx.x][threadIdx.y+i]; } }
}
__global__ void transposeCoarseGrained(float *odata, float *idata, int width, int height, int nreps)
{ __shared__ float block[TILE_DIM][TILE_DIM+1]; 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) { block[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] = block[threadIdx.y+i][threadIdx.x]; } }
}请注意，细粒度和粗粒度内核并不是实际的转置，因为在这两种情况下，odata都不是数据的转置，但正如您将看到的，它们在分析性能瓶颈时很有用。我们将这两种情况的性能结果添加到下表中:
细粒度转置具有与共享内存副本相似的性能，而粗粒度转置具有与合并转置和银行无冲突转置大致相同的性能。因此，性能瓶颈在于将数据写入全局内存中的转置位置。正如共享内存性能会因银行冲突而降低一样，通过分区露营进行全局内存访问也会导致类似的性能降低，这是我们接下来要研究的。
分区冲突
正如共享内存被划分为16个32位宽度的组一样，全局内存被划分为6个256字节宽度的分区(在8系列和9系列gpu上)或8个256字节宽度的分区(在200系列和10系列gpu上)。我们之前讨论过，为了有效地使用共享内存，半warp内的线程应该访问不同的银行，以便这些访问可以同时发生。如果半warp内的线程仅通过几个bank访问共享内存，则会发生bank冲突。
为了有效地使用全局内存，所有活动warp对全局内存的并发访问应该在各个分区之间平均分配。术语分区冲突用于描述全局内存访问定向通过分区子集的情况，导致请求在某些分区排队，而其他分区未使用。
合并关注的是半翘曲中的全局内存访问，而分区冲突关注的是活动半翘曲中的全局内存访问。由于分区冲突涉及活动线程块的行为，因此如何在多处理器上调度线程块的问题很重要。当内核启动时，块分配给多处理器的顺序由一维块ID决定，定义为:
bid = blockIdx.x + gridDim.x*blockIdx.y;这是网格中块的行-主排序。一旦达到最大占用率，就会根据需要将额外的块分配给多处理器。块完成的速度和顺序无法确定，因此活动块最初是连续的，但随着内核执行的进展，它们变得不那么连续。
如果我们回到矩阵转置并查看2048 x2048矩阵中的块如何映射到GTX 280上的分区，如下图所示，我们立即发现分区冲突是一个问题。
对于8个256字节宽度的分区，所有2048字节(或512个浮点数)的数据都映射到同一个分区。任何具有512列整数倍的浮点矩阵，例如我们的2048x2048矩阵，将包含其元素映射到单个分区的列。对于32 x32浮点数(或128 x 128字节)的块(其一维块id如图所示)，块的前两列中的所有数据都映射到同一个分区，对于其他对的块列也是如此(假设矩阵与分区段对齐)。
结合矩阵元素映射到分区的方式，以及块的调度方式，我们可以看到并发块将按行访问数据中的块，这些数据将大致均匀地分布在分区中，然而这些块将按列访问odata中的块，而odata通常只通过几个分区访问全局内存。
在将这个问题诊断为分区冲突之后，现在的问题是可以对此做些什么。与共享内存一样填充也是一个选项。向odata添加额外的64列(一个分区宽度)将导致一个tile的行依次映射到不同的分区。然而，对于某些应用程序来说，这种填充可能会变得令人望而却步。有一种更简单的解决方案，本质上涉及重新调度块的执行方式。
对角块重排
虽然程序员不能直接控制调度块的顺序(这是由自动内核变量blockldx的值决定的)，但程序员在如何解释blockldx的组件方面确实具有灵活性。给定组件blockldx的命名方式，
也就是x和y，人们通常认为这些分量指的是笛卡尔坐标系。然而，这并不一定是事实，人们可以选择其他方式。在笛卡尔解释中，可以交换这两个组件的角色，这将消除写入odata时的分区露营问题，但是这只是将问题转移到从数据中读取数据。
在读取数据和写入odata时避免分区露营的一种方法是对blockldx的组件使用对角线解释:它们组件表示通过矩阵的瓷砖的不同对角线切片，x组件表示沿着每个对角线的距离。对于4x4块矩阵，在下图的顶部显示了blockldx组件的笛卡尔和对角线解释，以及在底部产生的一维块ID。
在我们讨论在矩阵转置中使用blockldx分量的对角解释的优点之前，我们简要地提到如何使用坐标映射有效地实现它。这种技术在编写新内核时很有用，但在修改现有内核以使用对角线(或其他)对blockldx字段的解释时更是如此。如果blockldx。x和blockldx。Y表示对角线坐标，则(对于块方阵)对应的笛卡尔坐标由以下映射给出:
blockIdx_y = blockIdx.x;
blockIdx_x = (blockIdx.x+blockIdx.y)%gridDim.x;只需在内核的开头包含前两行代码，并假设对blockldx字段进行笛卡尔解释来编写内核，只是使用blockldx_x和blockldx_y代替blockldx。x和blockldx。Y，在整个内核中。这正是下面的转置对角线内核所做的:
__global__ void transposeDiagonal(float *odata, float *idata, int width, int height, int nreps)
{ __shared__ float tile[TILE_DIM][TILE_DIM+1]; int blockIdx_x, blockIdx_y; // diagonal reordering if (width == height) { blockIdx_y = blockIdx.x; blockIdx_x = (blockIdx.x+blockIdx.y)%gridDim.x;} else { int bid = blockIdx.x + gridDim.x*blockIdx.y; blockIdx_y = bid%gridDim.y; blockIdx_x = ((bid/gridDim.y)+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]; } }
}这里我们允许使用方阵和非方阵。一般情况下可以使用非方阵的映射，但是，方阵的简单表达式计算速度更快，在适当的时候更可取。
如果我们重新查看下图中的2048x2048矩阵，我们可以看到对角线重新排序是如何解决分区露营问题的。在对角线的情况下，当从数据读取数据并向odata写入数据时，就像从数据读取数据时在笛卡尔的情况下一样，成对的磁片循环通过分区。
下表中对角核的性能反映了这一点。:在内核内对全局内存的读写进行循环时测量的带宽当在内核上循环时，性能会略有下降，这可能是由在共享内存副本的几个百分点之内。于在计算blockldx x和
blockldx_y。然而，即使有这种性能下降，对角转置的带宽是其他完全转置的四倍以上。
Summary
在本文中，我们通过一系列逐步优化的转置内核讨论了GPU内存管理的几个方面。该序列是使用CUDA进行性能调优的典型序列。提高有效带宽的第一步是确保全局内存访问是合并的，这可以将性能提高一个数量级。
第二步是查看共享内存库冲突。在本研究中，消除共享内存库冲突似乎对性能几乎没有影响，但这主要是由于它与其他优化相关的应用:存储库冲突的影响被分区冲突掩盖了。通过在对角线重新排序的转置中删除共享内存数组的填充，可以看到银行冲突对性能有相当大的影响。
虽然合并和库冲突将随着问题大小的变化而保持相对一致，但分区冲突取决于问题大小并且在不同的硬件世代中有所不同。本例中特定大小的矩阵在基于g80的卡上由于分区数量不同(8系列上有6个分区，而200系列上有8个分区)而导致的分区冲突导致的性能下降要小得多。
转置内核的最终版本绝不代表可以实现的最高水平的优化。Tile大小、每个线程的元素数量和指令优化都可以提高转置内核和复制内核的性能。但在这项研究中，我们只关注了影响最大的问题。
Appendix A - Host Code
#include <stdio.h>
// kernels transpose/copy a tile of TILE_DIM x TILE_DIM elements
// using a TILE_DIM x BLOCK_ROWS thread block, so that each thread
// transposes TILE_DIM/BLOCK_ROWS elements. TILE_DIM must be an
// integral multiple of BLOCK_ROWS
#define TILE_DIM 32
#define BLOCK_ROWS 8
// Number of repetitions used for timing.
#define NUM_REPS 100
int
main( int argc, char** argv)
{ // set matrix size const int size_x = 2048, size_y = 2048; // kernel pointer and descriptor void (*kernel)(float *, float *, int, int, int); char *kernelName; // execution configuration parameters dim3 grid(size_x/TILE_DIM, size_y/TILE_DIM), threads(TILE_DIM,BLOCK_ROWS); // CUDA events cudaEvent_t start, stop; // size of memory required to store the matrix const int mem_size = sizeof(float) * size_x*size_y; // allocate host memory float *h_idata = (float*) malloc(mem_size); float *h_odata = (float*) malloc(mem_size); float *transposeGold = (float *) malloc(mem_size); float *gold; // allocate device memory float *d_idata, *d_odata; cudaMalloc( (void**) &d_idata, mem_size); cudaMalloc( (void**) &d_odata, mem_size); // initalize host data for(int i = 0; i < (size_x*size_y); ++i) h_idata[i] = (float) i; // copy host data to device cudaMemcpy(d_idata, h_idata, mem_size, cudaMemcpyHostToDevice );
// Compute reference transpose solution computeTransposeGold(transposeGold, h_idata, size_x, size_y); // print out common data for all kernels printf("\nMatrix size: %dx%d, tile: %dx%d, block: %dx%d\n\n", size_x, size_y, TILE_DIM, TILE_DIM, TILE_DIM, BLOCK_ROWS); printf("Kernel\t\t\tLoop over kernel\tLoop within kernel\n"); printf("------\t\t\t----------------\t------------------\n"); // // loop over different kernels // for (int k = 0; k<8; k++) { // set kernel pointer switch (k) { case 0: kernel = © kernelName = "simple copy "; break;case 1: kernel = ©SharedMem; kernelName = "shared memory copy "; break;case 2: kernel = &transposeNaive; kernelName = "naive transpose "; break;case 3: kernel = &transposeCoalesced; kernelName = "coalesced transpose "; break;case 4: kernel = &transposeNoBankConflicts; kernelName = "no bank conflict trans"; break;case 5: kernel = &transposeCoarseGrained; kernelName = "coarse-grained "; break;case 6: kernel = &transposeFineGrained; kernelName = "fine-grained "; break;case 7: kernel = &transposeDiagonal; kernelName = "diagonal transpose "; break;} // set reference solution // NB: fine- and coarse-grained kernels are not full // transposes, so bypass check if (kernel == © || kernel == ©SharedMem) { gold = h_idata; } else if (kernel == &transposeCoarseGrained || kernel == &transposeFineGrained) { gold = h_odata; } else { gold = transposeGold; } // initialize events, EC parameters cudaEventCreate(&start); cudaEventCreate(&stop); // warmup to avoid timing startup
kernel<<<grid, threads>>>(d_odata, d_idata, size_x,size_y, 1); // take measurements for loop over kernel launches cudaEventRecord(start, 0); for (int i=0; i < NUM_REPS; i++) { kernel<<<grid, threads>>>(d_odata, d_idata,size_x,size_y,1); } cudaEventRecord(stop, 0); cudaEventSynchronize(stop); float outerTime; cudaEventElapsedTime(&outerTime, start, stop); cudaMemcpy(h_odata,d_odata, mem_size, cudaMemcpyDeviceToHost); int res = comparef(gold, h_odata, size_x*size_y); if (res != 1) printf("*** %s kernel FAILED ***\n", kernelName); // take measurements for loop inside kernel cudaEventRecord(start, 0); kernel<<<grid,threads>>> (d_odata, d_idata, size_x, size_y, NUM_REPS); cudaEventRecord(stop, 0); cudaEventSynchronize(stop); float innerTime; cudaEventElapsedTime(&innerTime, start, stop); cudaMemcpy(h_odata,d_odata, mem_size, cudaMemcpyDeviceToHost); res = comparef(gold, h_odata, size_x*size_y); if (res != 1) printf("*** %s kernel FAILED ***\n", kernelName); // report effective bandwidths float outerBandwidth = 2.*1000*mem_size/(1024*1024*1024)/(outerTime/NUM_REPS); float innerBandwidth = 2.*1000*mem_size/(1024*1024*1024)/(innerTime/NUM_REPS); printf("%s\t%5.2f GB/s\t\t%5.2f GB/s\n", kernelName, outerBandwidth, innerBandwidth);} // cleanup free(h_idata); free(h_odata); free(transposeGold); cudaFree(d_idata); cudaFree(d_odata); cudaEventDestroy(start); cudaEventDestroy(stop); return 0;
}
