在CUDA中优化矩阵转置

Polaris北极星少女

已于 2024-08-30 19:33:41 修改

阅读量1.1k

点赞数 18

分类专栏： GPU编程文章标签：机器学习算法人工智能

于 2024-08-20 18:51:32 首次发布

本文链接：https://blog.csdn.net/qq_62704693/article/details/141358466

版权

GPU编程专栏收录该内容

7 篇文章 1 订阅

订阅专栏

https://dmacssite.github.io/materials/MatrixTranspose.pdfhttps://dmacssite.github.io/materials/MatrixTranspose.pdf

Chapter 1. Introduction

矩阵转置优化CUDA内存管理

本文档讨论了CUDA应用程序性能的各个方面，这些方面与有效使用GPU内存和应用于矩阵转置的数据管理有关。特别地，本文档讨论了以下内存使用问题:

合并数据传输到和从全局内存
共享内存bank冲突
分区冲突

主机和设备之间的数据传输，以及常量和纹理存储器。这里没有讨论高效内存使用的其他方面，例如合并和分区冲突都处理全局设备和片上内存之间的数据传输，而共享内存库冲突处理片上共享内存。
读者应该熟悉基本的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);

// 记录开始事件，作为测量的起始时间
cudaEventRecord(start, 0); 

// 循环执行内核函数 NUM_REPS 次
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);

// 记录开始事件，作为测量的起始时间
cudaEventRecord(start, 0); 

// 启动内核函数一次，并传递 NUM_REPS 作为循环次数参数
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__`开头是为了表明这是一个在 GPU 上执行的全局函数，也称为内核函数。以下是其具体含义： **一、执行位置** - 通常情况下，普通的 C/C++函数在 CPU 上执行。而以`__global__`开头的函数则是在 GPU 上执行的特殊函数。当调用这样的函数时，实际上是启动了一组 GPU 线程来并行执行这个函数。 **二、调用方式** - 这种内核函数不能像普通函数那样直接调用，而是需要通过特定的方式从主机（CPU）端调用，例如使用 CUDA 的运行时 API（如`cudaLaunchKernel`）或者使用 C++的模板函数调用语法（如`<<<gridSize, blockSize>>>`）。 **三、并行执行** - GPU 具有大量的并行处理单元，当启动一个内核函数时，多个线程会同时执行这个函数。每个线程都有自己独立的执行路径，并且可以通过线程索引（如这里的`threadIdx.x`、`threadIdx.y`、`blockIdx.x`、`blockIdx.y`）来访问不同的数据和执行不同的操作。 **四、参数传递** - 内核函数可以接收参数，这些参数通常是指向数据在 GPU 内存中的指针、维度信息（如这里的`width`、`height`）以及其他控制参数（如`nreps`）。这些参数从主机端传递到 GPU 端，以便内核函数能够正确地处理数据。总的来说，以`__global__`开头的函数是 CUDA 编程中用于在 GPU 上进行并行计算的关键组成部分，它允许开发者利用 GPU 的强大并行处理能力来加速计算密集型任务。

__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]; 
 } 
 } 
}

__global__ void copy(float *odata, float* idata, int width, int height, int nreps) 
{ 
    // 计算当前线程在x轴上的索引
    int xIndex = blockIdx.x * TILE_DIM + threadIdx.x; 
    
    // 计算当前线程在y轴上的索引
    int yIndex = blockIdx.y * TILE_DIM + threadIdx.y; 
    
    // 计算当前线程处理的数据在一维数组中的索引位置
    int index = xIndex + width * yIndex; 
    
    // 外层循环，重复nreps次
    for (int r = 0; r < nreps; r++) { 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行处理
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将idata数组中的数据拷贝到odata数组中
            odata[index + i * width] = idata[index + i * width]; 
        } 
    } 
}

#include <cuda_runtime.h>
#include <iostream>

#define TILE_DIM 32       // 每个线程块的维度
#define BLOCK_ROWS 32      // 每个线程块中的行数

// 核函数定义
__global__ void copy(float *odata, float* idata, int width, int height, int nreps) 
{ 
    // 计算当前线程在x轴上的索引
    int xIndex = blockIdx.x * TILE_DIM + threadIdx.x; 
    
    // 计算当前线程在y轴上的索引
    int yIndex = blockIdx.y * TILE_DIM + threadIdx.y; 
    
    // 计算当前线程处理的数据在一维数组中的索引位置
    int index = xIndex + width * yIndex; 
    
    // 将idata数组中的数据拷贝到odata数组中
    if (xIndex < width && yIndex < height) {
        odata[index] = idata[index];
    }
}

// 主机代码
int main() {
    int width = 1024;    // 矩阵的宽度
    int height = 1024;   // 矩阵的高度
    int nreps = 100;     // 重复次数

    // 分配主机内存
    float *h_idata = new float[width * height];
    float *h_odata = new float[width * height];

    // 初始化输入数据
    for(int i = 0; i < width * height; i++) {
        h_idata[i] = static_cast<float>(i);
    }

    // 分配设备内存
    float *d_idata, *d_odata;
    cudaMalloc(&d_idata, width * height * sizeof(float));
    cudaMalloc(&d_odata, width * height * sizeof(float));

    // 拷贝输入数据到设备
    cudaMemcpy(d_idata, h_idata, width * height * sizeof(float), cudaMemcpyHostToDevice);

    // 设置CUDA网格和块的维度
    dim3 dimBlock(TILE_DIM, BLOCK_ROWS);
    dim3 dimGrid((width + TILE_DIM - 1) / TILE_DIM, (height + TILE_DIM - 1) / TILE_DIM);

    // 调用核函数
    copy<<<dimGrid, dimBlock>>>(d_odata, d_idata, width, height, nreps);

    // 将结果从设备拷贝回主机
    cudaMemcpy(h_odata, d_odata, width * height * sizeof(float), cudaMemcpyDeviceToHost);

    // 检查结果（简化的检查方式）
    bool success = true;
    for(int i = 0; i < width * height; i++) {
        if(h_idata[i] != h_odata[i]) {
            success = false;
            break;
        }
    }

    if (success) {
        std::cout << "Data copy was successful!" << std::endl;
    } else {
        std::cout << "Data copy failed!" << std::endl;
    }

    // 释放内存
    delete[] h_idata;
    delete[] h_odata;
    cudaFree(d_idata);
    cudaFree(d_odata);

    return 0;
}

在这段代码中，`#define BLOCK_ROWS 8`定义了一个预处理常量。其含义如下： **一、在核函数中的作用** 1. 在核函数`copy`中，这个常量用于控制内层循环的步长以及确定每个线程块在处理数据时的行为。具体来说，在内层循环`for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS)`中，它决定了每次循环迭代时处理数据的间隔。例如，这里每次循环会将数据的索引增加`BLOCK_ROWS`的值，从而以固定的步长进行数据处理。 2. 这个常量也影响了核函数中数据的拷贝操作。通过以特定的步长进行处理，可以更有效地利用线程块中的线程来并行处理数据，提高数据拷贝的效率。 **二、整体性能优化方面的意义** 1. 通过调整`BLOCK_ROWS`的值，可以优化 GPU 上的线程块的组织和数据处理方式。不同的 GPU 架构和问题规模可能需要不同的`BLOCK_ROWS`值来达到最佳性能。较小的值可能会导致更多的线程启动，但每个线程处理的数据量较少；较大的值可能会减少线程启动的开销，但可能会导致某些线程闲置，因为数据量可能不足以充分利用所有线程。 2. 合理设置`BLOCK_ROWS`可以更好地平衡线程块的大小和数据处理的粒度，以充分利用 GPU 的并行计算能力，提高程序的整体性能。

这两种时间的区别在于内核启动的开销，这在不同的内核之间应该是一致的，以及在每个内核开始时计算矩阵索引所花费的时间。此外，在内核启动上循环也可以作为一种同步机制。当内核从宿主代码中的循环中多次启动时，一个内核启动的所有块必须在下一次启动的任何块开始之前完成执行。因此，每次循环选代都会重置活动块集和内存访问模式。当在内核内执行循环时，活动线程块集在执行定时循环过程中有更多的机会分散。
计时代码的两种方法都提供了有用的度量方法，第一种方法指示通常使用什么作为总体性能度量，第二种方法用于比较内核之间的数据移动时间。
在下一节中，我们将介绍从主机代码调用的不同内核，每个内核都解决不同的性能问题。本研究中的所有内核都启动尺寸为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上的循环用于多次计时数据从输入到输出数组的传输。

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];
 } 
 } 
}

__global__ void transposeNaive(float *odata, float* idata, int width, int height, int nreps) 
{ 
    // 计算当前线程在x轴上的索引
    int xIndex = blockIdx.x * TILE_DIM + threadIdx.x; 
    
    // 计算当前线程在y轴上的索引
    int yIndex = blockIdx.y * TILE_DIM + threadIdx.y; 
    
    // 计算输入数据数组中当前线程处理的数据的位置索引
    int index_in = xIndex + width * yIndex; 
    
    // 计算输出数据数组中当前线程处理的数据的位置索引
    int index_out = yIndex + height * xIndex; 
    
    // 外层循环，重复nreps次
    for (int r = 0; r < nreps; r++) { 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行处理
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将输入数组idata中的数据拷贝到输出数组odata中，进行转置
            odata[index_out + i] = idata[index_in + i * width];
        } 
    } 
}

#include <cuda_runtime.h>
#include <iostream>

#define TILE_DIM 32       // 每个线程块的维度
#define BLOCK_ROWS 8      // 每个线程块中的行数

// CUDA kernel 函数：实现矩阵转置
__global__ void transposeNaive(float *odata, float* idata, int width, int height, int nreps) 
{ 
    // 计算当前线程在 x 轴上的索引
    int xIndex = blockIdx.x * TILE_DIM + threadIdx.x; 
    
    // 计算当前线程在 y 轴上的索引
    int yIndex = blockIdx.y * TILE_DIM + threadIdx.y; 
    
    // 计算输入数组中当前线程处理的数据的位置索引
    int index_in = xIndex + width * yIndex; 
    
    // 计算输出数组中当前线程处理的数据的位置索引
    int index_out = yIndex + height * xIndex; 
    
    // 外层循环，重复 nreps 次
    for (int r = 0; r < nreps; r++) { 
        
        // 内层循环，通过步长 BLOCK_ROWS 对数据进行处理
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将输入数组 idata 中的数据拷贝到输出数组 odata 中，进行转置
            odata[index_out + i] = idata[index_in + i * width];
        } 
    } 
}

// 主机代码
int main() {
    int width = 1024;    // 矩阵的宽度
    int height = 1024;   // 矩阵的高度
    int nreps = 100;     // 重复次数

    // 分配主机内存
    float *h_idata = new float[width * height];
    float *h_odata = new float[width * height];

    // 初始化输入数据
    for(int i = 0; i < width * height; i++) {
        h_idata[i] = static_cast<float>(i);
    }

    // 分配设备内存
    float *d_idata, *d_odata;
    cudaMalloc(&d_idata, width * height * sizeof(float));
    cudaMalloc(&d_odata, width * height * sizeof(float));

    // 拷贝输入数据到设备
    cudaMemcpy(d_idata, h_idata, width * height * sizeof(float), cudaMemcpyHostToDevice);

    // 设置CUDA网格和块的维度
    dim3 dimBlock(TILE_DIM, BLOCK_ROWS);
    dim3 dimGrid((width + TILE_DIM - 1) / TILE_DIM, (height + TILE_DIM - 1) / TILE_DIM);

    // 调用核函数
    transposeNaive<<<dimGrid, dimBlock>>>(d_odata, d_idata, width, height, nreps);

    // 将结果从设备拷贝回主机
    cudaMemcpy(h_odata, d_odata, width * height * sizeof(float), cudaMemcpyDeviceToHost);

    bool isCorrect = true;

    // 遍历输出数据并验证是否与预期结果匹配
    for (int x = 0; x < width; ++x) {
        for (int y = 0; y < height; ++y) {
            int inputIndex = x + width * y;
            int outputIndex = y + height * x;
            if (h_odata[outputIndex] != h_idata[inputIndex]) {
                isCorrect = false;
                std::cerr << "数据不匹配！位置: (" << x << ", " << y << ") 输入值: " << h_idata[inputIndex]
                          << " 输出值: " << h_odata[outputIndex] << std::endl;
            }
        }
    }

    if (isCorrect) {
        std::cout << "转置结果验证成功！" << std::endl;
    } else {
        std::cout << "转置结果验证失败！" << std::endl;
    }

    // 释放内存
    delete[] h_idata;
    delete[] h_odata;
    cudaFree(d_idata);
    cudaFree(d_odata);

    return 0;
}

几乎与上面的复制内核相同，只是index(用于访问复制内核的输入和输出数组中的元素的数组索引)被两个索引index in(相当于复制内核中的index)和index_out取代。每个执行内核的线程将四个元素从输入矩阵的一列转置到它们在输出矩阵的一行中的转置位置。

这两个内核在2048 x2048矩阵上使用GTX280的性能如下表所示:

#include <iostream>
#include <cuda_runtime.h>

#define TILE_DIM 16
#define BLOCK_ROWS 16

// CUDA内核：简单复制
__global__ void copy(float *odata, const 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 + yIndex * width;

    for (int r = 0; r < nreps; r++) {
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) {
            if (xIndex < width && yIndex < height) {
                odata[index + i * width] = idata[index + i * width];
            }
        }
    }
}

// CUDA内核：朴素矩阵转置
__global__ void transposeNaive(float *odata, const 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) {
            if (xIndex < width && yIndex < height) {
                odata[index_out + i * height] = idata[index_in + i * width];
            }
        }
    }
}

// 主机代码：初始化矩阵，调用内核，验证结果
int main() {
    int width = 32;
    int height = 32;
    int nreps = 10;

    // 分配内存
    size_t size = width * height * sizeof(float);
    float *h_idata = new float[width * height];
    float *h_odata = new float[width * height];
    float *d_idata, *d_odata;

    cudaMalloc(&d_idata, size);
    cudaMalloc(&d_odata, size);

    // 初始化输入矩阵
    for (int i = 0; i < width * height; i++) {
        h_idata[i] = static_cast<float>(i);
    }

    // 拷贝数据到设备
    cudaMemcpy(d_idata, h_idata, size, cudaMemcpyHostToDevice);

    // 设置线程块和网格大小
    dim3 threadsPerBlock(TILE_DIM, TILE_DIM);
    dim3 numBlocks((width + TILE_DIM - 1) / TILE_DIM, (height + TILE_DIM - 1) / TILE_DIM);

    // 调用简单复制内核
    copy<<<numBlocks, threadsPerBlock>>>(d_odata, d_idata, width, height, nreps);

    // 拷贝结果回主机
    cudaMemcpy(h_odata, d_odata, size, cudaMemcpyDeviceToHost);

    // 验证复制结果
    bool correct = true;
    for (int i = 0; i < width * height; i++) {
        if (h_odata[i] != h_idata[i]) {
            correct = false;
            break;
        }
    }
    std::cout << "Copy result is " << (correct ? "correct" : "incorrect") << std::endl;

    // 打印复制结果矩阵的一部分
    std::cout << "Copy result (partial):" << std::endl;
    for (int y = 0; y < std::min(height, 10); y++) {
        for (int x = 0; x < std::min(width, 10); x++) {
            std::cout << h_odata[y * width + x] << " ";
        }
        std::cout << std::endl;
    }

    // 调用朴素转置内核
    transposeNaive<<<numBlocks, threadsPerBlock>>>(d_odata, d_idata, width, height, nreps);

    // 拷贝结果回主机
    cudaMemcpy(h_odata, d_odata, size, cudaMemcpyDeviceToHost);

    // 验证转置结果
    correct = true;
    for (int y = 0; y < height; y++) {
        for (int x = 0; x < width; x++) {
            if (h_odata[x * height + y] != h_idata[y * width + x]) {
                correct = false;
                break;
            }
        }
        if (!correct) break;
    }
    std::cout << "Transpose result is " << (correct ? "correct" : "incorrect") << std::endl;

    // 打印转置结果矩阵的一部分
    std::cout << "Transpose result (partial):" << std::endl;
    for (int y = 0; y < std::min(height, 10); y++) {
        for (int x = 0; x < std::min(width, 10); x++) {
            std::cout << h_odata[y * width + x] << " ";
        }
        std::cout << std::endl;
    }

    // 释放内存
    delete[] h_idata;
    delete[] h_odata;
    cudaFree(d_idata);
    cudaFree(d_odata);

    return 0;
}


### 1. **CUDA 内核函数定义**

代码定义了两个CUDA内核函数：`copy` 和 `transposeNaive`。我们重点关注 `transposeNaive` 内核函数，它负责执行矩阵的朴素转置操作。

### 2. **主机代码初始化**

在主机代码（`main()` 函数）中，以下步骤初始化并准备执行CUDA内核：

- **矩阵参数设置**：
  - 矩阵的宽度 `width` 和高度 `height` 被设置为 32。
  - `nreps` 参数用于指定内核在性能测量中的重复执行次数。

- **内存分配**：
  - 主机端分配两个大小为 `width * height` 的一维浮点数组 `h_idata` 和 `h_odata`。
  - 设备端（GPU）分配与主机相同大小的两个数组 `d_idata` 和 `d_odata`。

- **矩阵初始化**：
  - `h_idata` 数组被初始化为从0到(宽度*高度-1)的线性序列。这代表输入矩阵的每个元素。

- **数据传输**：
  - 使用 `cudaMemcpy` 函数将主机端 `h_idata` 的数据拷贝到设备端 `d_idata`，以便GPU内核函数使用。

### 3. **设置线程块和网格**

- **线程块和网格维度**：
  - 线程块的维度 `dim3 threadsPerBlock(TILE_DIM, TILE_DIM)` 定义了每个块的线程数量。这里 `TILE_DIM` 被设置为 16，因此每个线程块是一个 16x16 的网格。
  - 网格维度 `dim3 numBlocks((width + TILE_DIM - 1) / TILE_DIM, (height + TILE_DIM - 1) / TILE_DIM)` 定义了总的线程块数量。这里 `width` 和 `height` 都是 32，因此网格的每个维度都是 2。

### 4. **CUDA 内核函数执行**

- **调用朴素转置内核函数**：
  - 使用 `<<<numBlocks, threadsPerBlock>>>` 调用 `transposeNaive` 内核函数。每个线程块包含 256 个线程（16x16），整个网格包含 4 个线程块（2x2）。

- **CUDA 内核执行细节**：

  1. **线程索引计算**：
     - 每个线程都有一个唯一的线程索引，它通过 `blockIdx.x`, `blockIdx.y`, `threadIdx.x`, 和 `threadIdx.y` 计算得到。
     - `xIndex` 和 `yIndex` 分别计算了当前线程在整个网格中的X轴和Y轴的全局索引。

  2. **计算输入和输出矩阵的索引**：
     - `index_in = xIndex + width * yIndex` 计算了输入矩阵中当前线程负责处理的元素的位置索引。
     - `index_out = yIndex + height * xIndex` 计算了输出矩阵中转置后的位置索引。

  3. **数据处理循环**：
     - 内层循环 `for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS)` 允许线程块中的每个线程在Y轴方向上处理多行数据，而不仅仅是一行。
     - `odata[index_out + i * height] = idata[index_in + i * width]` 将输入矩阵中的元素拷贝到输出矩阵中相应的转置位置。
     - 外层循环 `for (int r = 0; r < nreps; r++)` 用于重复执行转置操作 `nreps` 次，以便测量性能。

  4. **边界检查**：
     - 使用 `if (xIndex < width && yIndex + i < height)` 确保线程的索引不会超出矩阵的边界。这样可以避免访问无效内存区域，确保程序安全运行。

### 5. **结果回传与验证**

- **结果传回主机**：
  - 内核执行完毕后，使用 `cudaMemcpy` 函数将结果矩阵 `d_odata` 从设备端拷贝回主机端的 `h_odata`。

- **结果验证**：
  - 通过检查 `h_odata` 中的数据与预期的转置结果进行对比，验证转置操作是否正确完成。

- **输出部分结果**：
  - 为了方便检查，代码打印了部分矩阵的内容。这部分代码是为了调试方便，可以帮助确认内核操作的正确性。

### 6. **资源释放**

- **释放内存**：
  - 主机和设备端的内存被释放，确保程序不会造成内存泄漏。

### 7. **代码执行流程总结**

整个过程可以概括为：

1. **初始化**：主机分配内存并初始化数据。
2. **传输数据**：将输入数据从主机传输到设备。
3. **内核执行**：使用CUDA内核函数进行矩阵的朴素转置操作。
4. **传回结果**：将结果从设备传回主机。
5. **验证结果**：验证转置操作是否正确。
6. **输出结果**：打印部分矩阵内容以确认正确性。
7. **释放资源**：释放主机和设备的内存。

这个过程演示了一个CUDA程序从启动到结束的完整执行流程。朴素转置方法尽管简单，但可以作为理解CUDA编程的基础，有助于进一步优化和复杂的CUDA应用程序的开发。

#include <iostream>
#include <cuda_runtime.h>

#define TILE_DIM 16
#define BLOCK_ROWS 8
#define NUM_REPS 1000

__global__ void transposeNaive(float *odata, const 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) {
            if (xIndex < width && yIndex + i < height) {
                odata[index_out + i * height] = idata[index_in + i * width];
            }
        }
    }
}

int main() {
    int size_x = 32;
    int size_y = 32;
    int mem_size = size_x * size_y * sizeof(float);

    float *h_idata = (float *)malloc(mem_size);
    float *h_odata = (float *)malloc(mem_size);

    for (int i = 0; i < size_x * size_y; i++) {
        h_idata[i] = (float)i;
    }

    float *d_idata, *d_odata;
    cudaMalloc((void**)&d_idata, mem_size);
    cudaMalloc((void**)&d_odata, mem_size);
    cudaMemcpy(d_idata, h_idata, mem_size, cudaMemcpyHostToDevice);

    dim3 threadsPerBlock(TILE_DIM, TILE_DIM);
    dim3 numBlocks((size_x + TILE_DIM - 1) / TILE_DIM, (size_y + TILE_DIM - 1) / TILE_DIM);

    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);

    // 内核外部循环测量
    cudaEventRecord(start, 0);
    for (int i = 0; i < NUM_REPS; i++) {
        transposeNaive<<<numBlocks, threadsPerBlock>>>(d_odata, d_idata, size_x, size_y, 1);
    }
    cudaEventRecord(stop, 0);
    cudaEventSynchronize(stop);
    float outerTime;
    cudaEventElapsedTime(&outerTime, start, stop);

    // 内核内部循环测量
    cudaEventRecord(start, 0);
    transposeNaive<<<numBlocks, threadsPerBlock>>>(d_odata, d_idata, size_x, size_y, NUM_REPS);
    cudaEventRecord(stop, 0);
    cudaEventSynchronize(stop);
    float innerTime;
    cudaEventElapsedTime(&innerTime, start, stop);

    // 打印时间测量结果
    std::cout << "Outer loop time: " << outerTime << " ms" << std::endl;
    std::cout << "Inner loop time: " << innerTime << " ms" << std::endl;

    cudaMemcpy(h_odata, d_odata, mem_size, cudaMemcpyDeviceToHost);

    // 打印部分结果用于验证
    for (int i = 0; i < size_x * size_y; i++) {
        std::cout << h_odata[i] << " ";
        if ((i + 1) % size_x == 0) std::cout << std::endl;
    }

    // 释放资源
    cudaFree(d_idata);
    cudaFree(d_odata);
    free(h_idata);
    free(h_odata);

    cudaEventDestroy(start);
    cudaEventDestroy(stop);

    return 0;
}

复制核和原始转置核之间代码的微小差异对性能有深远的影响--几乎有两个数量级的影响。这就引出了我们的第一个优化技术:全局内存合并。

合并转置

由于设备内存比片上内存具有更高的延迟和更低的带宽，因此必须特别注意如何执行全局内存访问，在我们的示例中，从data加载数据并将数据存储在odata中。如果满足某些条件，由half-warp线程进行的所有全局内存访问可以合并到一个或两个事务中。这些标准取决于设备的计算能力，这可以通过运行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个内存事务。

在CUDA编程中，**warp**（译为“波前”）是一个非常重要的概念。它指的是CUDA设备上的一组并行执行的线程。具体来说，一个warp通常由32个线程组成，这些线程会同时执行相同的指令集。

### 详细解释
- **同步执行**：在warp内的所有线程会同步执行相同的指令，也就是说，这32个线程在同一时刻执行同一个操作，但它们可以处理不同的数据。
- **SIMT架构**：CUDA采用了SIMT（Single Instruction, Multiple Threads）架构，也就是单一指令、多线程。warp是这种架构的基本执行单位，因此CUDA会一次处理一个warp的线程。
- **线程分组**：在编写CUDA程序时，通常会将线程分成块（blocks），每个块可以包含多个warp。warp的大小是固定的32个线程，这也是CUDA硬件的基础特性。

### 性能相关
- **分支和发散**：如果warp内的线程需要执行不同的指令，比如在条件语句中发生了分支（有些线程需要执行“if”部分，而有些需要执行“else”部分），warp内的线程会出现“发散”（divergence），这会降低执行效率，因为不同的指令需要依次执行，而不能并行。
- **内存访问模式**：warp内的线程在访问全局内存时，如果能做到合并访问（coalesced access），即这些线程一起访问连续的内存地址，会极大地提高内存访问效率。如果访问不连续，称为未合并访问（uncoalesced access），则会导致性能下降。

总的来说，warp是CUDA并行计算的基本单元，理解warp的工作原理对优化CUDA程序性能非常重要。

简而言之，如果内存访问合并到计算能力为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]; 
 } 
 } 
}

__global__ void transposeCoalesced(float *odata, float *idata, int width, int height, int nreps) 
{ 
    // 在共享内存中定义一个TILE_DIM x TILE_DIM的二维数组tile，用于存储矩阵块
    __shared__ float tile[TILE_DIM][TILE_DIM]; 
    
    // 计算输入数据在x轴和y轴的索引
    int xIndex = blockIdx.x * TILE_DIM + threadIdx.x; 
    int yIndex = blockIdx.y * TILE_DIM + threadIdx.y; 
    
    // 计算输入数据数组中当前线程处理的数据的位置索引
    int index_in = xIndex + (yIndex) * width; 
    
    // 重新计算xIndex和yIndex，以适应输出数据的转置
    xIndex = blockIdx.y * TILE_DIM + threadIdx.x; 
    yIndex = blockIdx.x * TILE_DIM + threadIdx.y; 
    
    // 计算输出数据数组中当前线程处理的数据的位置索引
    int index_out = xIndex + (yIndex) * height; 
    
    // 外层循环，重复nreps次
    for (int r = 0; r < nreps; r++) { 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行处理
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将输入数组idata中的数据加载到共享内存tile中
            tile[threadIdx.y + i][threadIdx.x] = idata[index_in + i * width]; 
        } 
        
        // 等待所有线程完成共享内存的数据加载
        __syncthreads(); 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行转置和写回
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将共享内存tile中的数据写入输出数组odata
            odata[index_out + i * height] = tile[threadIdx.x][threadIdx.y + i]; 
        } 
    } 
}

#include <iostream>
#include <cstdio>
#include <cuda_runtime.h>

#define TILE_DIM 16        // 定义每个块的维度
#define BLOCK_ROWS 4       // 定义每块中的行数

// CUDA kernel 函数：实现矩阵转置
__global__ void transposeCoalesced(float *odata, float *idata, int width, int height, int nreps) 
{ 
    // 在共享内存中定义一个TILE_DIM x TILE_DIM的二维数组tile，用于存储矩阵块
    __shared__ float tile[TILE_DIM][TILE_DIM]; 
    
    // 计算输入数据在x轴和y轴的索引
    int xIndex = blockIdx.x * TILE_DIM + threadIdx.x; 
    int yIndex = blockIdx.y * TILE_DIM + threadIdx.y; 
    
    // 计算输入数据数组中当前线程处理的数据的位置索引
    int index_in = xIndex + (yIndex) * width; 
    
    // 重新计算xIndex和yIndex，以适应输出数据的转置
    xIndex = blockIdx.y * TILE_DIM + threadIdx.x; 
    yIndex = blockIdx.x * TILE_DIM + threadIdx.y; 
    
    // 计算输出数据数组中当前线程处理的数据的位置索引
    int index_out = xIndex + (yIndex) * height; 
    
    // 外层循环，重复nreps次
    for (int r = 0; r < nreps; r++) { 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行处理
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将输入数组idata中的数据加载到共享内存tile中
            tile[threadIdx.y + i][threadIdx.x] = idata[index_in + i * width]; 
        } 
        
        // 等待所有线程完成共享内存的数据加载
        __syncthreads(); 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行转置和写回
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将共享内存tile中的数据写入输出数组odata
            odata[index_out + i * height] = tile[threadIdx.x][threadIdx.y + i]; 
        } 
    } 
}


// 检查 CUDA 错误的函数
void checkCudaErrors(cudaError_t err) {
    if (err != cudaSuccess) {
        std::cerr << "CUDA 错误: " << cudaGetErrorString(err) << std::endl;
        exit(EXIT_FAILURE);
    }
}

// 验证转置结果的函数
void verifyTranspose(float *hostInput, float *hostOutput, int width, int height) {
    bool isCorrect = true;

    // 遍历输出数据并验证是否与预期结果匹配
    for (int x = 0; x < width; ++x) {
        for (int y = 0; y < height; ++y) {
            int inputIndex = x + width * y;
            int outputIndex = y + height * x;
            if (hostOutput[outputIndex] != hostInput[inputIndex]) {
                isCorrect = false;
                std::cerr << "数据不匹配！位置: (" << x << ", " << y << ") 输入值: " << hostInput[inputIndex]
                          << " 输出值: " << hostOutput[outputIndex] << std::endl;
            }
        }
    }

    if (isCorrect) {
        std::cout << "转置结果验证成功！" << std::endl;
    } else {
        std::cout << "转置结果验证失败！" << std::endl;
    }
}

int main() {
    const int width = 256;     // 矩阵的宽度
    const int height = 256;    // 矩阵的高度
    const int size = width * height * sizeof(float); // 数据大小
    const int nreps = 10;      // 循环重复次数

    float *h_idata, *h_odata;  // 主机端数据指针
    float *d_idata, *d_odata;  // 设备端数据指针

    // 分配主机内存
    h_idata = (float*)malloc(size);
    h_odata = (float*)malloc(size);

    // 初始化主机数据
    for (int i = 0; i < width * height; ++i) {
        h_idata[i] = static_cast<float>(i);
    }

    // 分配设备内存
    checkCudaErrors(cudaMalloc((void**)&d_idata, size));
    checkCudaErrors(cudaMalloc((void**)&d_odata, size));

    // 将数据从主机拷贝到设备
    checkCudaErrors(cudaMemcpy(d_idata, h_idata, size, cudaMemcpyHostToDevice));

    // 定义网格和块的维度
    dim3 threads(TILE_DIM, BLOCK_ROWS);     // 每个块中的线程数
    dim3 grid((width + TILE_DIM - 1) / TILE_DIM, (height + TILE_DIM - 1) / TILE_DIM); // 网格中的块数

    // 启动 kernel 函数
    transposeCoalesced<<<grid, threads>>>(d_odata, d_idata, width, height, nreps);

    // 等待所有的 CUDA 操作完成
    checkCudaErrors(cudaDeviceSynchronize());

    // 将结果从设备拷贝回主机
    checkCudaErrors(cudaMemcpy(h_odata, d_odata, size, cudaMemcpyDeviceToHost));

    // 验证结果
    verifyTranspose(h_idata, h_odata, width, height);

    // 释放设备内存
    checkCudaErrors(cudaFree(d_idata));
    checkCudaErrors(cudaFree(d_odata));

    // 释放主机内存
    free(h_idata);
    free(h_odata);

    return 0;
}

简单复制和朴素转置操作中的内存合并情况

1. **简单复制**：对于简单复制操作，在上述讨论的任何计算能力的设备上，从`idata`的所有加载操作都能合并。在`i`循环的每次迭代中，每个线程束读取 16 个连续的 32 位字，即一个瓦片（tile）的半行。通过`cudaMalloc()`分配设备内存并选择`TILE_DIM`为 16 的倍数可以确保与内存段对齐，因此所有加载都是合并的。在向`odata`写入时，在`i`循环的每次迭代中，一个线程束以合并的方式写入一个tile'的半行。

2. **朴素转置**：在朴素转置操作中，向`odata`写入时情况不同。在`i`循环的每次迭代中，一个线程束将一半的列浮点数写入不同的内存段，无论计算能力如何，都会导致 16 个单独的内存事务。

避免非合并全局内存访问的方法

避免非合并全局内存访问的方法是将数据读入共享内存，然后让每个线程束在共享内存中访问非连续的位置，以便将连续的数据写入`odata`。在共享内存中进行非连续访问没有像在全局内存中那样的性能损失。但是，这个过程要求一个瓦片（tile）中的每个元素由不同的线程访问，因此需要一个`__synchthreads()`调用，以确保在从共享内存写入`odata`之前，所有从`idata`到共享内存的读取都已完成。

下面给出了合并转置内核中半个线程束（half warp）的数据流描述。半个线程束将`idata`矩阵瓦片的四行半数据写入由黄色线段指示的共享内存中 32x32 的数组“tile”。在调用`__syncthreads()`以确保对“tile”的所有写入都已完成后，半个线程束将“tile”中的四列半数据写入`odata`矩阵瓦片的四行半中，由绿色线段指示。

在CUDA编程中，描述了“合并转置（coalesced transpose）”内核中半个warp的数据流。

### 数据流动过程：
1. **写入共享内存**：
- 半个warp负责将`idata`矩阵tile的四个半行数据写入共享内存中的32x32数组“tile”。这些写入操作由黄色线段表示。由于warp的合并访问特性，这些数据的写入是高效且合并的。

2. **同步线程**：
- 写入共享内存后，会调用`__syncthreads()`函数，以确保所有线程的写入操作都已经完成。这一步非常重要，因为共享内存中的数据需要在所有线程都完成写入后才能被正确读取和操作。

3. **从共享内存写入输出数据**：
- 在同步完成后，半个warp从共享内存的“tile”数组中读取四个半列数据，并将它们写入`odata`矩阵tile的四个半行。这个操作由绿色线段表示。

这个过程有效地解决了矩阵转置中可能出现的未合并内存访问问题，通过利用共享内存来重新排列数据，从而在最终写入全局内存时实现了合并访问，提升了执行效率。

为了详细说明“合并转置”（coalesced transpose）内核中半个warp的数据流，我们可以通过一个具体的示例来讲解这个过程。

### 假设场景：

- **warp**：一个warp包含32个线程，因此半个warp包含16个线程。
- **`idata`矩阵**：假设`idata`是一个64x64大小的矩阵，我们只关注其中的一个32x32的tile。
- **共享内存tile数组**：共享内存中分配了一个32x32大小的数组，称为“tile”。

### 数据流示例：

1. **写入共享内存：**
   - 假设我们处理的是`idata`矩阵的第一个tile，它在`idata`中的位置为`idata[0:31, 0:31]`。
   - 假设`idata`矩阵的前四行如下：

   ```
   idata:
   行0: 1 2 3 ... 16 ... 32
   行1: 33 34 35 ... 48 ... 64
   行2: 65 66 67 ... 80 ... 96
   行3: 97 98 99 ... 112 ... 128
   ```

   - 半个warp（16个线程）负责将这四行的前16个元素（即四个半行）写入共享内存中的tile数组的前四行。写入过程如下：
   
   ```
   tile数组（共享内存）:
   [0, 0] = idata[0, 0] = 1      [0, 1] = idata[0, 1] = 2     ... [0, 15] = idata[0, 15] = 16
   [1, 0] = idata[1, 0] = 33     [1, 1] = idata[1, 1] = 34    ... [1, 15] = idata[1, 15] = 48
   [2, 0] = idata[2, 0] = 65     [2, 1] = idata[2, 1] = 66    ... [2, 15] = idata[2, 15] = 80
   [3, 0] = idata[3, 0] = 97     [3, 1] = idata[3, 1] = 98    ... [3, 15] = idata[3, 15] = 112
   ```

   - 在这个过程中，由于每个线程负责处理连续的内存地址，数据的写入是合并访问的。

2. **同步线程：**
   - 调用`__syncthreads()`确保所有线程完成了对共享内存的写入。这样可以保证接下来所有的线程都能从共享内存中读取到完整、正确的数据。

3. **从共享内存写入输出数据：**
   - 同步完成后，半个warp从共享内存的tile数组中读取数据，并将其写入`odata`矩阵的相应位置。
   - 现在，半个warp将`tile`数组的前四列数据写入`odata`矩阵tile的前四行。这一操作可以描述如下：

   ```
   tile数组（共享内存）:
   [0, 0] = 1      [1, 0] = 33     [2, 0] = 65     [3, 0] = 97
   [0, 1] = 2      [1, 1] = 34     [2, 1] = 66     [3, 1] = 98
   ...
   [0, 15] = 16    [1, 15] = 48    [2, 15] = 80    [3, 15] = 112
   ```

   - 写入到`odata`的过程：

   ```
   odata:
   行0: 1 33 65 97 ... 
   行1: 2 34 66 98 ... 
   行2: 3 35 67 99 ... 
   行3: 4 36 68 100 ...
   ```

   - 这个写入过程会使得输出矩阵`odata`中的数据与`idata`的转置版本一致。
   - 由于在共享内存中数据已经重新排列，写入`odata`的过程可以通过合并访问完成，从而大大提高了内存访问效率。

### **总结：**
- 通过这个示例，可以看到“合并转置”如何利用共享内存来重新排列数据，使得最终写入全局内存的操作是合并访问的。这样既避免了未合并访问带来的性能损失，也实现了高效的矩阵转置。

通过改进odata对内存的访问模式，写操作被合并，我们看到了性能的提高

尽管“合并转置”相比“普通转置”（naïve transpose）在有效带宽上有显著提升，但其性能仍与“简单拷贝”存在较大的差距。分析表明，这种性能差距并非由于转置操作中额外的索引计算导致，因为在“内核循环”列（Loop in kernel）中的结果显示，即使将索引计算分摊到100次数据移动中，仍然存在较大的性能差异。

一个可能导致这种性能差距的原因是“合并转置”中所需的同步屏障（`__syncthreads()`）。为了验证这一点，可以通过以下利用共享内存并包含`__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]; 
 } 
 } 
}

__global__ void copySharedMem(float *odata, float *idata, int width, int height, int nreps) {
    // 在共享内存中定义一个TILE_DIM x TILE_DIM的二维数组tile
    __shared__ float tile[TILE_DIM][TILE_DIM]; 
    
    // 计算当前线程在x轴上的全局索引
    int xIndex = blockIdx.x * TILE_DIM + threadIdx.x; 
    
    // 计算当前线程在y轴上的全局索引
    int yIndex = blockIdx.y * TILE_DIM + threadIdx.y; 
    
    // 计算输入数据数组中当前线程处理的数据的位置索引
    int index = xIndex + width * yIndex; 
    
    // 外层循环，重复nreps次
    for (int r = 0; r < nreps; r++) { 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行处理
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将输入数组idata中的数据加载到共享内存tile中
            tile[threadIdx.y + i][threadIdx.x] = idata[index + i * width]; 
        } 
        
        // 等待所有线程完成共享内存的数据加载
        __syncthreads(); 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行处理
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将共享内存tile中的数据拷贝回输出数组odata中
            odata[index + i * width] = tile[threadIdx.y + i][threadIdx.x]; 
        } 
    } 
}

在这个内核中，__syncthreads()调用并不是为了确保内核成功执行而必须的，因为线程之间并不共享数据。它仅仅是为了评估合并转置中的同步屏障的开销而被包含在内。修改后的结果显示在下表中。

#include <cuda_runtime.h>
#include <iostream>

#define TILE_DIM 32       // 每个线程块的维度
#define BLOCK_ROWS 8      // 每个线程块中的行数

// 核函数定义
__global__ void copySharedMem(float *odata, float *idata, int width, int height, int nreps) {
    // 在共享内存中定义一个TILE_DIM x TILE_DIM的二维数组tile
    __shared__ float tile[TILE_DIM][TILE_DIM]; 
    
    // 计算当前线程在x轴上的全局索引
    int xIndex = blockIdx.x * TILE_DIM + threadIdx.x; 
    
    // 计算当前线程在y轴上的全局索引
    int yIndex = blockIdx.y * TILE_DIM + threadIdx.y; 
    
    // 计算输入数据数组中当前线程处理的数据的位置索引
    int index = xIndex + width * yIndex; 
    
    // 外层循环，重复nreps次
    for (int r = 0; r < nreps; r++) { 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行处理
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将输入数组idata中的数据加载到共享内存tile中
            tile[threadIdx.y + i][threadIdx.x] = idata[index + i * width]; 
        } 
        
        // 等待所有线程完成共享内存的数据加载
        __syncthreads(); 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行处理
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将共享内存tile中的数据拷贝回输出数组odata中
            odata[index + i * width] = tile[threadIdx.y + i][threadIdx.x]; 
        } 
    } 
}

// 主机代码
int main() {
    int width = 1024;    // 矩阵的宽度
    int height = 1024;   // 矩阵的高度
    int nreps = 100;     // 重复次数

    // 分配主机内存
    float *h_idata = new float[width * height];
    float *h_odata = new float[width * height];

    // 初始化输入数据
    for(int i = 0; i < width * height; i++) {
        h_idata[i] = static_cast<float>(i);
    }

    // 分配设备内存
    float *d_idata, *d_odata;
    cudaMalloc(&d_idata, width * height * sizeof(float));
    cudaMalloc(&d_odata, width * height * sizeof(float));

    // 拷贝输入数据到设备
    cudaMemcpy(d_idata, h_idata, width * height * sizeof(float), cudaMemcpyHostToDevice);

    // 设置CUDA网格和块的维度
    dim3 dimBlock(TILE_DIM, BLOCK_ROWS);
    dim3 dimGrid((width + TILE_DIM - 1) / TILE_DIM, (height + TILE_DIM - 1) / TILE_DIM);

    // 调用核函数
    copySharedMem<<<dimGrid, dimBlock>>>(d_odata, d_idata, width, height, nreps);

    // 将结果从设备拷贝回主机
    cudaMemcpy(h_odata, d_odata, width * height * sizeof(float), cudaMemcpyDeviceToHost);

    // 检查结果（简化的检查方式）
    bool success = true;
    for(int i = 0; i < width * height; i++) {
        if(h_idata[i] != h_odata[i]) {
            success = false;
            break;
        }
    }

    if (success) {
        std::cout << "Data copy was successful!" << std::endl;
    } else {
        std::cout << "Data copy failed!" << std::endl;
    }

    // 释放内存
    delete[] h_idata;
    delete[] h_odata;
    cudaFree(d_idata);
    cudaFree(d_odata);

    return 0;
}

共享内存复制的结果似乎表明，使用共享内存和同步屏障对性能影响很小，至少从“内核循环”（“Loop in kernel”）列中比较简单复制和共享内存复制时可以看出这一点。然而，当比较合并转置和共享内存复制内核时，有一个需要解决的性能瓶颈：共享内存的访问方式可能导致共享内存bank冲突

共享内存bank冲突

共享内存被分成16个大小相等的内存模块，称为存储库，这些存储库的组织方式是将连续的32位字分配给连续的存储库。这些bank可以同时被访问，为了获得最大的带宽进出共享内存，半曲的线程应该访问与不同bank相关联的共享内存。此规则的例外情况是，当半warp中的所有线程读取相同的共享内存地址时，这会导致广播，其中该地址的数据在一个事务中发送给半warp的所有线程。
在分析CUDA应用程序时，可以使用warpserialize标志来确定共享内存库冲突是否发生在任何内核中。一般来说，这个标志也反映了原子和常量内存的使用，但是在我们的示例中这两者都不存在。
合并转置使用32x32的浮点数共享内存数组。对于这个大小的数组，列k和k+16中的所有数据都映射到同一个库。因此，当从共享内存中的tile写入部分列到odata中的行时半warp会经历16路bank冲突并序列化请求。避免这种冲突的一个简单方法是将共享内存数组填充一列:

1. **共享内存bank**：
- 共享内存被分成16个大小相等的内存模块，称为（banks），这些bank以这样的方式组织：连续的32位字（words）被分配到连续的bank。
- 这些bank可以同时被访问，为了实现对共享内存的最大带宽，半个warp中的线程应该访问与不同bank相关联的共享内存。

2. **广播例外**：
- 如果半个warp中的所有线程读取相同的共享内存地址，就会发生广播（broadcast），该地址的数据会在一次事务中发送到半个warp的所有线程。

3. **检测bank冲突**：
- 在分析CUDA应用程序时，可以使用`warp_serialize`标志来确定内核中是否发生了共享内存bank冲突。虽然这个标志也反映了原子操作和常量内存的使用，但在这个示例中并不存在这些情况。

4. **合并转置中的bank冲突**：
- 合并转置使用了一个32x32大小的浮点数共享内存数组。对于这种大小的数组，列k和k+16的数据被映射到相同的bank。
- 结果是，当将共享内存中的部分列数据写入`odata`中的行时，半个warp会遇到16路bank冲突，导致请求的串行化。

5. **解决bank冲突的方法**：
- 通过在共享内存数组中添加一列的填充来避免这种冲突：
```cpp
__shared__ float tile[TILE_DIM][TILE_DIM+1];
```
- 这个填充不会影响半个warp写入共享内存时的bank访问模式，仍然是无冲突的，但现在当半个warp访问列中的数据时，也不会发生冲突。

6. **性能提升**：
- 通过添加一列的填充，内核的性能提升了，现在这个内核既是合并的，也是无bank冲突的，并且这一改进的性能已经添加到表格中。

总结：通过在共享内存数组中添加一列的填充，可以有效避免共享内存bank冲突，从而提升CUDA内核的性能。

__shared__ float tile[TILE_DIM][TILE_DIM+1];

添加的填充列并不会影响半个warp向共享内存写入数据时的共享内存银行访问模式，这种写入操作仍然是无冲突的。但是，通过添加一列填充后，当半个warp访问共享内存中的列数据时，也不会再发生冲突。现在这个内核既是合并访问的，也是无共享内存银行冲突的，其性能已经被添加到我们下面的表格中。

#include <iostream>
#include <cstdio>
#include <cuda_runtime.h>

#define TILE_DIM 16        // 定义每个块的维度
#define BLOCK_ROWS 4       // 定义每块中的行数

// CUDA kernel 函数：实现矩阵转置
__global__ void transposeNoBankConflicts(float *odata, float *idata, int width, int height, int nreps) 
{ 
    // 在共享内存中定义一个TILE_DIM x TILE_DIM的二维数组tile，用于存储矩阵块
    __shared__ float tile[TILE_DIM][TILE_DIM+1];
    
    // 计算输入数据在x轴和y轴的索引
    int xIndex = blockIdx.x * TILE_DIM + threadIdx.x; 
    int yIndex = blockIdx.y * TILE_DIM + threadIdx.y; 
    
    // 计算输入数据数组中当前线程处理的数据的位置索引
    int index_in = xIndex + (yIndex) * width; 
    
    // 重新计算xIndex和yIndex，以适应输出数据的转置
    xIndex = blockIdx.y * TILE_DIM + threadIdx.x; 
    yIndex = blockIdx.x * TILE_DIM + threadIdx.y; 
    
    // 计算输出数据数组中当前线程处理的数据的位置索引
    int index_out = xIndex + (yIndex) * height; 
    
    // 外层循环，重复nreps次
    for (int r = 0; r < nreps; r++) { 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行处理
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将输入数组idata中的数据加载到共享内存tile中
            tile[threadIdx.y + i][threadIdx.x] = idata[index_in + i * width]; 
        } 
        
        // 等待所有线程完成共享内存的数据加载
        __syncthreads(); 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行转置和写回
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将共享内存tile中的数据写入输出数组odata
            odata[index_out + i * height] = tile[threadIdx.x][threadIdx.y + i]; 
        } 
    } 
}
// 主机代码
int main() {
    int width = 1024;    // 矩阵的宽度
    int height = 1024;   // 矩阵的高度
    int nreps = 100;     // 重复次数

    // 分配主机内存
    float *h_idata = new float[width * height];
    float *h_odata = new float[width * height];

    // 初始化输入数据
    for(int i = 0; i < width * height; i++) {
        h_idata[i] = static_cast<float>(i);
    }

    // 分配设备内存
    float *d_idata, *d_odata;
    cudaMalloc(&d_idata, width * height * sizeof(float));
    cudaMalloc(&d_odata, width * height * sizeof(float));

    // 拷贝输入数据到设备
    cudaMemcpy(d_idata, h_idata, width * height * sizeof(float), cudaMemcpyHostToDevice);

    // 设置CUDA网格和块的维度
    dim3 dimBlock(TILE_DIM, BLOCK_ROWS);
    dim3 dimGrid((width + TILE_DIM - 1) / TILE_DIM, (height + TILE_DIM - 1) / TILE_DIM);

    // 调用核函数
    transposeNoBankConflicts<<<dimGrid, dimBlock>>>(d_odata, d_idata, width, height, nreps);

    // 将结果从设备拷贝回主机
    cudaMemcpy(h_odata, d_odata, width * height * sizeof(float), cudaMemcpyDeviceToHost);

    bool isCorrect = true;

    // 遍历输出数据并验证是否与预期结果匹配
    for (int x = 0; x < width; ++x) {
        for (int y = 0; y < height; ++y) {
            int inputIndex = x + width * y;
            int outputIndex = y + height * x;
            if (h_odata[outputIndex] != h_idata[inputIndex]) {
                isCorrect = false;
                std::cerr << "数据不匹配！位置: (" << x << ", " << y << ") 输入值: " << h_idata[inputIndex]
                          << " 输出值: " << h_odata[outputIndex] << std::endl;
            }
        }
    }

    if (isCorrect) {
        std::cout << "转置结果验证成功！" << std::endl;
    } else {
        std::cout << "转置结果验证失败！" << std::endl;
    }

    // 释放内存
    delete[] h_idata;
    delete[] h_odata;
    cudaFree(d_idata);
    cudaFree(d_odata);

    return 0;
}

### 1. 什么是共享内存银行冲突？

共享内存分为多个“银行”（通常为16个）。每个银行可以在一次时钟周期内处理一个请求，这意味着如果一个warp中的所有线程都访问不同的银行，所有请求可以同时完成。否则，如果多个线程同时访问同一个银行中的不同地址，就会发生“银行冲突”，这些请求必须被顺序处理，导致性能下降。

### 2. 数据示例说明银行冲突

假设我们有一个32x32的共享内存数组`tile[TILE_DIM][TILE_DIM]`，其中`TILE_DIM = 32`。每个元素是32位（即4字节）的浮点数（`float`类型）。共享内存中的元素分布如下：

- 第1列（即k列）中的数据映射到共享内存的第1个银行。
- 第17列（即k+16列）中的数据也映射到同一个银行。

```plaintext
银行0: tile[0][0], tile[0][16], tile[1][0], tile[1][16], ...
银行1: tile[0][1], tile[0][17], tile[1][1], tile[1][17], ...
...
银行15: tile[0][15], tile[0][31], tile[1][15], tile[1][31], ...
```

### 3. 发生银行冲突的情况

假设一个半warp（即16个线程）尝试从共享内存中读取第k列的数据，并将其写入全局内存的一个行中，同时另一个半warp读取k+16列的数据，并将其写入到全局内存的另一行中。因为第k列和第k+16列的数据都映射到相同的银行中，这样就会产生银行冲突：

```plaintext
银行0: tile[0][0], tile[0][16], tile[1][0], tile[1][16], ...
```

例如，假设第一个线程读取`tile[0][0]`，第二个线程读取`tile[0][16]`，由于这两个数据映射到同一个银行0，所以线程必须按顺序进行操作，而不是并行执行。这就会导致16路银行冲突。

### 4. 解决银行冲突的方法

通过在共享内存数组中增加一列的填充，可以打破这种冲突：

```cpp
__shared__ float tile[TILE_DIM][TILE_DIM+1];
```

增加的这列填充数据会使得连续的列数据不再映射到同一个银行。例如：

```plaintext
银行0: tile[0][0], tile[0][17], tile[1][0], tile[1][17], ...
```

由于现在第k列的数据与第k+16列的数据不再映射到同一个银行，线程可以并行访问这些数据，不会发生冲突。

### 5. 结果

通过在共享内存中添加额外的一列填充（`TILE_DIM+1`），可以消除银行冲突，从而提升性能。例如，如果原来由于银行冲突导致每次访问需要16个时钟周期，现在消除冲突后，可能只需要1个时钟周期，这将大大提高内核的执行速度。

**总结：** 数据说明了在共享内存访问过程中，由于不同列映射到同一个银行而产生的银行冲突问题。通过简单的添加一列填充，可以有效避免冲突，显著提高程序的执行效率。

虽然在共享内存数组中添加填充确实消除了共享内存银行冲突（这一点通过使用CUDA分析器检查`warp_serialize`标志得到了确认），但在这个阶段实现填充对性能的提升影响不大。因此，合并访问且无共享内存银行冲突的转置操作与共享内存拷贝之间仍然存在较大的性能差距。在接下来的部分中，我们将对转置操作进行分解，以确定导致性能下降的原因。

分解的转置

在上表中，最佳优化的转置操作与共享内存拷贝之间的性能差距超过了四倍。这不仅体现在对内核启动次数的循环测量中，还体现在内核内部循环100次时的测量结果中，尽管额外的索引计算成本在这100次迭代中已被分摊。

为了进一步研究这一差距，我们重新审视了转置操作的数据流，并将其与拷贝操作的数据流进行比较，这两者都在下面图表的上半部分进行了说明。基本上，拷贝代码和转置操作之间有两个主要区别：在tile（块）中转置数据，以及将数据写入转置后的tile。我们可以通过实现两个内核，每个内核单独执行其中一个操作来隔离这两个组件之间的性能差异。

如图表的下半部分所示，细粒度转置内核在tile内部转置数据，但将tile写入拷贝操作通常会写入的位置。而粗粒度转置内核则将tile写入odata矩阵中的转置位置，但不会在tile内部转置数据。

这两个内核的源代码如下:

__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]; 
 } 
 } 
}

#include <cuda_runtime.h>
#include <iostream>

#define TILE_DIM 32
#define BLOCK_ROWS 8

__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) { 
            if (xIndex < width && yIndex + i < height) {
                block[threadIdx.y + i][threadIdx.x] = idata[index + i * width]; 
            }
        } 
        
        __syncthreads(); 
        
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            if (xIndex < width && yIndex + i < height) {
                odata[index + i * height] = block[threadIdx.x][threadIdx.y + i]; 
            }
        } 
    } 
}

void checkCUDAError(const char *msg) {
    cudaError_t err = cudaGetLastError();
    if (err != cudaSuccess) {
        std::cerr << "CUDA error: " << msg << " : " << cudaGetErrorString(err) << std::endl;
        exit(EXIT_FAILURE);
    }
}

int main() {
    int width = 1024;
    int height = 1024;
    int nreps = 10;

    size_t size = width * height * sizeof(float);
    float *h_idata = new float[width * height];
    float *h_odata = new float[width * height];

    // Initialize input data
    for (int i = 0; i < width * height; ++i) {
        h_idata[i] = static_cast<float>(i);
    }

    float *d_idata, *d_odata;
    cudaMalloc(&d_idata, size);
    cudaMalloc(&d_odata, size);

    cudaMemcpy(d_idata, h_idata, size, cudaMemcpyHostToDevice);

    dim3 threads(TILE_DIM, BLOCK_ROWS);
    dim3 blocks((width + TILE_DIM - 1) / TILE_DIM, (height + TILE_DIM - 1) / TILE_DIM);
    
    transposeFineGrained<<<blocks, threads>>>(d_odata, d_idata, width, height, nreps);
    checkCUDAError("Kernel execution failed");

    cudaMemcpy(h_odata, d_odata, size, cudaMemcpyDeviceToHost);
    checkCUDAError("Memory copy failed");

    // Verification of results (simple print for the first 10 elements)
    for (int i = 0; i < 10; ++i) {
        std::cout << h_odata[i] << " ";
        if (i % width == width - 1)
            std::cout << std::endl;
    }

    delete[] h_idata;
    delete[] h_odata;
    cudaFree(d_idata);
    cudaFree(d_odata);

    return 0;
}

#include <cuda_runtime.h>
#include <iostream>

#define TILE_DIM 32
#define BLOCK_ROWS 8

__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) { 
            if (xIndex < width && yIndex + i < height) {
                block[threadIdx.y + i][threadIdx.x] = idata[index_in + i * width]; 
            }
        } 
        
        __syncthreads(); 
        
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            if (xIndex < height && yIndex + i < width) {
                odata[index_out + i * height] = block[threadIdx.y + i][threadIdx.x]; 
            }
        } 
    } 
}

void checkCUDAError(const char *msg) {
    cudaError_t err = cudaGetLastError();
    if (err != cudaSuccess) {
        std::cerr << "CUDA error: " << msg << " : " << cudaGetErrorString(err) << std::endl;
        exit(EXIT_FAILURE);
    }
}

int main() {
    int width = 1024;
    int height = 1024;
    int nreps = 10;

    size_t size = width * height * sizeof(float);
    float *h_idata = new float[width * height];
    float *h_odata = new float[width * height];

    // Initialize input data
    for (int i = 0; i < width * height; ++i) {
        h_idata[i] = static_cast<float>(i);
    }

    float *d_idata, *d_odata;
    cudaMalloc(&d_idata, size);
    cudaMalloc(&d_odata, size);

    cudaMemcpy(d_idata, h_idata, size, cudaMemcpyHostToDevice);

    dim3 threads(TILE_DIM, BLOCK_ROWS);
    dim3 blocks((width + TILE_DIM - 1) / TILE_DIM, (height + TILE_DIM - 1) / TILE_DIM);
    
    transposeCoarseGrained<<<blocks, threads>>>(d_odata, d_idata, width, height, nreps);
    checkCUDAError("Kernel execution failed");

    cudaMemcpy(h_odata, d_odata, size, cudaMemcpyDeviceToHost);
    checkCUDAError("Memory copy failed");

    // Verification of results (simple print for the first 10 elements)
    for (int i = 0; i < 10; ++i) {
        std::cout << h_odata[i] << " ";
        if (i % width == width - 1)
            std::cout << std::endl;
    }

    delete[] h_idata;
    delete[] h_odata;
    cudaFree(d_idata);
    cudaFree(d_odata);

    return 0;
}

请注意，细粒度和粗粒度内核并不是实际的转置，因为在这两种情况下，odata都不是数据的转置，但正如您将看到的，它们在分析性能瓶颈时很有用。我们将这两种情况的性能结果添加到下表中:

细粒度转置的性能与共享内存拷贝相似，而粗粒度转置的性能则大致与合并且无银行冲突的转置相当。因此，性能瓶颈在于将数据写入全局内存中的转置位置。就像共享内存的性能会因银行冲突而下降一样，全局内存访问也可能因分区竞争（partition camping）而出现类似的性能下降问题，我们将对此进行进一步研究。

分区冲突

正如共享内存被分成16个32位宽的bank一样，全局内存在不同系列的GPU上被划分为不同数量的分区：8系列和9系列GPU有6个分区，而200系列和10系列GPU有8个分区，每个分区的宽度是256字节。我们之前讨论过，为了有效使用共享内存，半warp内的线程应该访问不同的bank，这样这些访问就可以同时发生。如果半warp内的线程只访问共享内存的几个bank，那么就会发生bank冲突。

有效使用全局内存的关键是，所有活跃warp对全局内存的同时访问应该均匀地分配到各个分区中。"分区冲突"（partition camping）这个术语描述了一种情况，即当全局内存访问被导向通过部分分区时，会导致一些分区上的请求排队等待，而其他分区则未被充分利用。

合并（coalescing）关注的是半warp内部的全局内存访问，而分区冲突则关注的是不同活跃半warp之间的全局内存访问。由于分区冲突涉及到活跃线程块在多处理器上的行为，因此线程块在多处理器上的调度方式变得尤为重要。当一个内核被启动时，块被分配到多处理器的顺序由一维块ID决定。

bid = blockIdx.x + gridDim.x*blockIdx.y;

这是网格中块的行主序排列。一旦达到最大占用率，就会根据需要将额外的块分配给多处理器。但是，块完成的速度和顺序是无法确定的，因此活跃块最初是连续的，但随着内核执行的进展，它们变得越来越不连续。

如果我们回到矩阵转置问题，看看2048x2048矩阵的tile如何在GTX 280上的分区映射，就像下面的图表所示，我们立刻就能看出分区冲突是一个问题。

对于8个256字节宽度的分区，所有2048字节(或512个浮点数)的数据都映射到同一个分区。任何具有512列整数倍的浮点矩阵，例如我们的2048x2048矩阵，将包含其元素映射到单个分区的列。对于32 x32浮点数(或128 x 128字节)的块(其一维块id如图所示)，块的前两列中的所有数据都映射到同一个分区，对于其他对的块列也是如此(假设矩阵与分区段对齐)。
结合矩阵元素映射到分区的方式，以及块的调度方式，我们可以看到并发块将按行访问数据中的块，这些数据将大致均匀地分布在分区中，然而这些块将按列访问odata中的块，而odata通常只通过几个分区访问全局内存。
在将这个问题诊断为分区冲突之后，现在的问题是可以对此做些什么。与共享内存一样填充也是一个选项。向odata添加额外的64列(一个分区宽度)将导致一个tile的行依次映射到不同的分区。然而，对于某些应用程序来说，这种填充可能会变得令人望而却步。有一种更简单的解决方案，本质上涉及重新调度块的执行方式。

**一、全局内存的分区结构** 就像共享内存被划分为 16 个 32 位宽度的存储体一样，全局内存在 8 系列和 9 系列 GPU 上被分为 6 个分区，在 200 系列和 10 系列 GPU 上被分为 8 个分区，每个分区的宽度为 256 字节。 **二、有效使用共享内存与全局内存的方式对比** 1. 共享内存：为了有效地使用共享内存，半个线程束中的线程应该访问不同的存储体，以便这些访问可以同时进行。如果半个线程束中的线程只通过少数几个存储体访问共享内存，就会发生存储体冲突。 2. 全局内存：为了有效地使用全局内存，所有活跃线程束对全局内存的并发访问应该在分区之间均匀分配。当全局内存访问集中在一部分分区上，导致请求在某些分区排队，而其他分区未被使用时，就会出现“分区占用”（partition camping）的情况。 **三、分区占用与线程块调度的关系** 1. 分区占用与活跃的半个线程束之间的全局内存访问有关，而共享内存的合并则与半个线程束内的全局内存访问有关。 2. 当启动一个内核时，线程块分配给多处理器的顺序是由一维的块 ID 决定的，即`bid = blockIdx.x + gridDim.x*blockIdx.y`，这是网格中线程块的行优先顺序。一旦达到最大占用率，根据需要将额外的线程块分配给多处理器。线程块完成的速度和顺序是不确定的，所以在开始时活跃的线程块是连续的，但随着内核的执行，它们会变得不那么连续。 3. 以矩阵转置为例，观察 2048x2048 的矩阵中的瓦片如何映射到 GTX 280 上的分区。在具有 8 个 256 字节宽度分区的情况下，步长为 2048 字节（或 512 个浮点数）的所有数据都映射到同一个分区。任何具有整数倍 512 列的浮点矩阵，如 2048x2048 的矩阵，其列中的元素将映射到单个分区。对于 32x32 个浮点数的瓦片（或 128x128 字节），其一维块 ID 如图所示，瓦片的前两列中的所有数据都映射到同一个分区，其他瓦片列对也是如此（假设矩阵是按照分区段对齐的）。 4. 结合矩阵元素如何映射到分区以及线程块是如何调度的，可以看到并发的线程块将在行方向上访问`idata`中的瓦片，这些访问将大致均匀地分布在各个分区中，但这些线程块将在列方向上访问`odata`中的瓦片，这通常会通过少数几个分区访问全局内存，从而导致分区占用问题。 **四、解决分区占用问题的方法** 1. 诊断出问题是分区占用后，一个解决方法是像处理共享内存一样进行填充。向`odata`添加额外的 64 列（一个分区的宽度）将使瓦片的行依次映射到不同的分区。然而，这种填充对于某些应用可能会变得不可行。 2. 有一个更简单的解决方案，主要涉及重新安排线程块的执行顺序。

为什么整数倍会映射到相同分区

当矩阵的列数是整数倍时（如 `2048` 列），为什么会导致这些列的数据映射到相同的内存分区呢？这个问题的核心在于全局内存的地址计算方式以及分区映射的机制。

### 1. 内存地址与分区映射
在 CUDA 中，GPU 的全局内存是按照一定的分区（partition）来管理的。每个分区的宽度通常是固定的（例如，256 字节）。当一个线程访问全局内存时，内存地址的低位决定了数据将映射到哪个分区。

假设全局内存被分为 8 个分区，每个分区的宽度为 256 字节。地址空间会根据分区宽度进行循环映射。具体来说：

- 地址 `0x00000000` 到 `0x000000FF` 会映射到分区 0。
- 地址 `0x00000100` 到 `0x000001FF` 会映射到分区 1。
- ...
- 地址 `0x00000700` 到 `0x000007FF` 会映射到分区 7。
- 地址 `0x00000800` 到 `0x000008FF` 又会回到分区 0。

这种循环映射的模式依赖于地址的低位字节数，而分区的数量（如 8 个）和每个分区的宽度（如 256 字节）决定了映射的周期性。

### 2. 为什么整数倍的步幅会导致相同映射
在矩阵运算中，步幅（stride）是指访问同一列的相邻元素在内存地址中的间隔。对于一个 `2048x2048` 的矩阵：

- 假设每个浮点数占 4 字节，则每行占用的内存大小为 `2048 × 4 字节 = 8192 字节`。
- 也就是说，相邻两行的同一列元素在内存中的地址相差 8192 字节。

由于 8192 字节是 256 字节（一个分区宽度）的 `32` 倍，因此 `2048` 列矩阵中相邻两行的同一列元素在内存中的地址低位是完全相同的。这导致相同列中的元素会不断地映射到同一个分区。例如：

- 第 1 列的第 1 个元素地址为 `0x00000000`，映射到分区 0。
- 第 1 列的第 2 个元素地址为 `0x00002000`，同样映射到分区 0（因为地址低位部分一致）。

这种情况发生在所有列上，所有元素都会重复映射到相同的分区。这就是为什么矩阵的列数是某个特定整数倍（例如 `2048` 列）时会导致所有列映射到相同分区的原因。

### 3. 如何避免这种情况
为了避免这种映射问题，我们可以：

- **改变步幅**：通过增加矩阵的列数，使得步幅不再是分区宽度的整数倍，从而打破这种周期性映射。
- **数据填充**：例如，在矩阵末尾添加额外的列，使得每行的内存占用不再是 `8192` 字节（即 `2048 × 4 字节`），从而使得内存地址的低位在每一行发生变化。

### 4. 结论
当矩阵的步幅是全局内存分区宽度的整数倍时，由于地址低位字节的一致性，所有列的元素会重复映射到相同的分区。这会导致分区阻塞，从而降低内存访问的效率。通过调整步幅或添加数据填充，可以有效避免这种情况。

对角块重排

虽然程序员无法直接控制线程块的调度顺序（该顺序由自动内核变量 `blockIdx` 的值决定），但程序员可以灵活地解释 `blockIdx` 的各个组成部分。鉴于 `blockIdx` 的分量命名为 `x` 和 `y`，通常人们会认为这些分量对应于笛卡尔坐标系。然而，这并非必须如此，程序员可以选择其他方式进行解释。在笛卡尔坐标系的解释中，可以交换这两个分量的角色，从而消除写入 `odata` 时的分区阻塞问题，但这只会将问题转移到读取 `idata` 上。

为了解决在读取 `idata` 和写入 `odata` 时的分区阻塞问题，可以采用对角线的解释方式：将 `blockIdx.y` 解释为矩阵中不同对角线切片的标识，而 `blockIdx.x` 则表示沿每条对角线的距离。图示展示了 `4x4` 块矩阵中，`blockIdx` 分量在笛卡尔坐标系和对角线解释中的分布方式，以及其结果的一维块 ID。

在讨论在矩阵转置中使用 `blockIdx` 分量的对角线解释的优点之前，我们先简要介绍如何通过坐标映射高效地实现这一解释。这种技术在编写新内核时很有用，但在修改现有内核以使用 `blockIdx` 字段的对角线（或其他）解释时更为重要。如果 `blockIdx.x` 和 `blockIdx.y` 表示对角线坐标，那么对于块矩阵来说，相应的笛卡尔坐标可以通过以下映射给出：

blockIdx_y = blockIdx.x;
blockIdx_x = (blockIdx.x + blockIdx.y) % gridDim.x;

只需在内核的开头加入上述两行代码，然后在整个内核中假设 `blockIdx` 字段的笛卡尔解释，并分别用 `blockIdx_x` 和 `blockIdx_y` 代替 `blockIdx.x` 和 `blockIdx.y`。这正是下面的 `transposeDiagonal` 内核中所做的操作。

__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]; 
 } 
 } 
}

__global__ void transposeDiagonal(float *odata, float *idata, int width, int height, int nreps) 
{ 
    // 在共享内存中定义一个TILE_DIM x (TILE_DIM+1)的二维数组tile
    __shared__ float tile[TILE_DIM][TILE_DIM+1]; 
    
    int blockIdx_x, blockIdx_y; 
    
    // 对角线重排序逻辑
    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; 
    } 
    
    // 计算当前线程在x轴上的全局索引
    int xIndex = blockIdx_x * TILE_DIM + threadIdx.x; 
    
    // 计算当前线程在y轴上的全局索引
    int yIndex = blockIdx_y * TILE_DIM + threadIdx.y; 
    
    // 计算输入数据数组中当前线程处理的数据的位置索引
    int index_in = xIndex + (yIndex) * width; 
    
    // 重新计算xIndex和yIndex，以适应输出数据的转置
    xIndex = blockIdx_y * TILE_DIM + threadIdx.x; 
    yIndex = blockIdx_x * TILE_DIM + threadIdx.y; 
    
    // 计算输出数据数组中当前线程处理的数据的位置索引
    int index_out = xIndex + (yIndex) * height; 
    
    // 外层循环，重复nreps次
    for (int r = 0; r < nreps; r++) { 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行处理
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将输入数组idata中的数据加载到共享内存tile中
            tile[threadIdx.y + i][threadIdx.x] = idata[index_in + i * width]; 
        } 
        
        // 等待所有线程完成共享内存的数据加载
        __syncthreads(); 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行处理
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将共享内存tile中的数据写回到输出数组odata中
            odata[index_out + i * height] = tile[threadIdx.x][threadIdx.y + i]; 
        } 
    } 
}

在这里，我们允许有方阵和非方阵。非方阵的映射可以在一般情况下使用，然而，方阵的更简单表达式计算速度更快，在适当的时候更可取。如果我们回顾下面图中的 2048x2048 矩阵，我们可以看到对角线重排序如何解决分区拥挤问题。在对角线情况下，从`idata`读取数据并写入`odata`时，就像在笛卡尔情况下从`idata`读取数据时一样，成对的分块在分区中循环。

好的，接下来我将详细解释每个步骤，包括各个线程块和线程的工作流程，如何在共享内存中存储数据，以及最终如何将数据转置并存储到输出矩阵 `odata` 中。

### 1. 矩阵和参数的设置

假设矩阵 `idata` 是一个 \(4 \times 4\) 的矩阵，数据如下：

```text
idata = [1, 2, 3, 4,
         5, 6, 7, 8,
         9, 10, 11, 12,
         13, 14, 15, 16]
```

`TILE_DIM = 2`，`BLOCK_ROWS = 1`，`nreps = 1`。

网格的大小为 \(2 \times 2\)，因此有 4 个线程块，每个线程块大小为 \(2 \times 2\)，每个线程块有 4 个线程。

### 2. 线程块分配及对角线重排序

因为这是一个方阵（宽度等于高度），所以我们使用方阵的对角线重排序逻辑。假设网格的原始排列为：

```text
(0,0)  (1,0)
(0,1)  (1,1)
```

经过对角线重排序后，网格中的线程块在 X 轴上的新索引为：

```text
(0,0)  (1,0)
(1,1)  (0,1)
```

- 线程块 (0,0) 保持不变，处理左上角的块。
- 线程块 (1,0) 保持不变，处理右上角的块。
- 线程块 (0,1) 和 (1,1) 交换，分别处理左下角和右下角的块。

### 3. 线程块内的工作

每个线程块处理一个 \(2 \times 2\) 的子块。以下详细解释每个线程块中的线程如何工作。

#### **线程块 (0, 0)**

- **线程 (0,0)**:
  - `xIndex = 0 * 2 + 0 = 0`
  - `yIndex = 0 * 2 + 0 = 0`
  - 从 `idata[0 + 0 * 4] = idata[0]` 读取数据，将 `1` 加载到 `tile[0][0]`。
  - 写回到 `odata[0 + 0 * 4] = odata[0]`，存储 `1`。

- **线程 (0,1)**:
  - `xIndex = 0 * 2 + 1 = 1`
  - `yIndex = 0 * 2 + 0 = 0`
  - 从 `idata[1 + 0 * 4] = idata[1]` 读取数据，将 `2` 加载到 `tile[0][1]`。
  - 写回到 `odata[4 + 0 * 4] = odata[4]`，存储 `5`。

- **线程 (1,0)**:
  - `xIndex = 0 * 2 + 0 = 0`
  - `yIndex = 1 * 2 + 0 = 2`
  - 从 `idata[0 + 2 * 4] = idata[8]` 读取数据，将 `9` 加载到 `tile[1][0]`。
  - 写回到 `odata[1 + 2 * 4] = odata[9]`，存储 `2`。

- **线程 (1,1)**:
  - `xIndex = 1 * 2 + 1 = 1`
  - `yIndex = 1 * 2 + 1 = 2`
  - 从 `idata[1 + 2 * 4] = idata[9]` 读取数据，将 `10` 加载到 `tile[1][1]`。
  - 写回到 `odata[5 + 1 * 4] = odata[10]`，存储 `6`。

共享内存 `tile` 的内容：

```text
tile = [1, 2,
        5, 6]
```

写回到 `odata` 中后部分输出为：

```text
odata = [1, 5, _, _,
         _, _, _, _,
         2, 6, _, _,
         _, _, _, _]
```

#### **线程块 (1, 0)**

- **线程 (0,0)**:
  - `xIndex = 1 * 2 + 0 = 2`
  - `yIndex = 0 * 2 + 0 = 0`
  - 从 `idata[2 + 0 * 4] = idata[2]` 读取数据，将 `3` 加载到 `tile[0][0]`。
  - 写回到 `odata[8 + 0 * 4] = odata[8]`，存储 `9`。

- **线程 (0,1)**:
  - `xIndex = 2 * 2 + 1 = 3`
  - `yIndex = 0 * 2 + 0 = 0`
  - 从 `idata[3 + 0 * 4] = idata[3]` 读取数据，将 `4` 加载到 `tile[0][1]`。
  - 写回到 `odata[12 + 0 * 4] = odata[12]`，存储 `13`。

- **线程 (1,0)**:
  - `xIndex = 1 * 2 + 0 = 2`
  - `yIndex = 1 * 2 + 1 = 2`
  - 从 `idata[2 + 2 * 4] = idata[10]` 读取数据，将 `11` 加载到 `tile[1][0]`。
  - 写回到 `odata[9 + 2 * 4] = odata[10]`，存储 `11`。

- **线程 (1,1)**:
  - `xIndex = 1 * 2 + 1 = 3`
  - `yIndex = 1 * 2 + 1 = 2`
  - 从 `idata[3 + 2 * 4] = idata[11]` 读取数据，将 `12` 加载到 `tile[1][1]`。
  - 写回到 `odata[13 + 3 * 4] = odata[15]`，存储 `14`。

共享内存 `tile` 的内容：

```text
tile = [3, 4,
        7, 8]
```

写回到 `odata` 中后部分输出为：

```text
odata = [1, 5, 9, 13,
         2, 6, 10, 14,
         3, 7, 11, 15,
         4, 8, 12, 16]
```

### 4. 最终结果

最终，所有线程块都完成了处理，`odata` 中存储的矩阵将是 `idata` 的转置：

```text
odata = [1, 5, 9, 13,
         2, 6, 10, 14,
         3, 7, 11, 15,
         4, 8, 12, 16]
```

这样，通过每个线程块内的共享内存 `tile`，数据被先读取到共享内存中，再通过同步操作后，将数据写回到全局内存 `odata`，并完成了矩阵的转置操作。

下面表格中的对角内核的性能反映了这一点。在内核中循环读取和写入全局内存时测量的带宽与共享内存副本的带宽相差几个百分点。当在核上循环时，性能略有下降，可能是由于计算`blockIdx_x`和`blockIdx_y`涉及到额外的计算。然而，即使有这种性能下降，对角线转置的带宽也是其他完全转置的四倍多。

#include <iostream>
#include <cstdio>
#include <cuda_runtime.h>

#define TILE_DIM 32       // 定义每个块的维度
#define BLOCK_ROWS 8       // 定义每块中的行数

// CUDA kernel 函数：实现矩阵转置
__global__ void transposeDiagonal(float *odata, float *idata, int width, int height, int nreps) 
{ 
    // 在共享内存中定义一个TILE_DIM x (TILE_DIM+1)的二维数组tile
    __shared__ float tile[TILE_DIM][TILE_DIM+1]; 
    
    int blockIdx_x, blockIdx_y; 
    
    // 对角线重排序逻辑
    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; 
    } 
    
    // 计算当前线程在x轴上的全局索引
    int xIndex = blockIdx_x * TILE_DIM + threadIdx.x; 
    
    // 计算当前线程在y轴上的全局索引
    int yIndex = blockIdx_y * TILE_DIM + threadIdx.y; 
    
    // 计算输入数据数组中当前线程处理的数据的位置索引
    int index_in = xIndex + (yIndex) * width; 
    
    // 重新计算xIndex和yIndex，以适应输出数据的转置
    xIndex = blockIdx_y * TILE_DIM + threadIdx.x; 
    yIndex = blockIdx_x * TILE_DIM + threadIdx.y; 
    
    // 计算输出数据数组中当前线程处理的数据的位置索引
    int index_out = xIndex + (yIndex) * height; 
    
    // 外层循环，重复nreps次
    for (int r = 0; r < nreps; r++) { 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行处理
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将输入数组idata中的数据加载到共享内存tile中
            tile[threadIdx.y + i][threadIdx.x] = idata[index_in + i * width]; 
        } 
        
        // 等待所有线程完成共享内存的数据加载
        __syncthreads(); 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行处理
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将共享内存tile中的数据写回到输出数组odata中
            odata[index_out + i * height] = tile[threadIdx.x][threadIdx.y + i]; 
        } 
    } 
}


// 主机代码
int main() {
    int width = 1024;    // 矩阵的宽度
    int height = 1024;   // 矩阵的高度
    int nreps = 100;     // 重复次数

    // 分配主机内存
    float *h_idata = new float[width * height];
    float *h_odata = new float[width * height];

    // 初始化输入数据
    for(int i = 0; i < width * height; i++) {
        h_idata[i] = static_cast<float>(i);
    }

    // 分配设备内存
    float *d_idata, *d_odata;
    cudaMalloc(&d_idata, width * height * sizeof(float));
    cudaMalloc(&d_odata, width * height * sizeof(float));

    // 拷贝输入数据到设备
    cudaMemcpy(d_idata, h_idata, width * height * sizeof(float), cudaMemcpyHostToDevice);

    // 设置CUDA网格和块的维度
    dim3 dimBlock(TILE_DIM, BLOCK_ROWS);
    dim3 dimGrid((width + TILE_DIM - 1) / TILE_DIM, (height + TILE_DIM - 1) / TILE_DIM);

    // 调用核函数
    transposeDiagonal<<<dimGrid, dimBlock>>>(d_odata, d_idata, width, height, nreps);

    // 将结果从设备拷贝回主机
    cudaMemcpy(h_odata, d_odata, width * height * sizeof(float), cudaMemcpyDeviceToHost);

    bool isCorrect = true;

    // 遍历输出数据并验证是否与预期结果匹配
    for (int x = 0; x < width; ++x) {
        for (int y = 0; y < height; ++y) {
            int inputIndex = x + width * y;
            int outputIndex = y + height * x;
            if (h_odata[outputIndex] != h_idata[inputIndex]) {
                isCorrect = false;
                std::cerr << "数据不匹配！位置: (" << x << ", " << y << ") 输入值: " << h_idata[inputIndex]
                          << " 输出值: " << h_odata[outputIndex] << std::endl;
            }
        }
    }

    if (isCorrect) {
        std::cout << "转置结果验证成功！" << std::endl;
    } else {
        std::cout << "转置结果验证失败！" << std::endl;
    }

    // 释放内存
    delete[] h_idata;
    delete[] h_odata;
    cudaFree(d_idata);
    cudaFree(d_odata);

    return 0;
}

Summary

在本文中，我们通过一系列逐步优化的转置内核讨论了 GPU 内存管理的几个方面。这个序列是使用 CUDA 进行性能调优的典型过程。提高有效带宽的第一步是确保全局内存访问是合并的，这可以将性能提高一个数量级。第二步是查看共享内存存储体冲突。在这项研究中，消除共享内存存储体冲突似乎对性能影响不大，然而，这在很大程度上是由于它在与其他优化相关的时间点上被应用：存储体冲突的影响被分区拥挤所掩盖。通过在对角线重排序的转置中去除共享内存数组的填充，可以看到存储体冲突对性能有相当大的影响。虽然随着问题规模的变化，合并和存储体冲突将保持相对一致，但分区拥挤取决于问题规模，并且在不同代的硬件上有所不同。在这个例子中的特定大小的矩阵，由于分区数量的不同，在基于 G80 的显卡上由于分区拥挤而导致的性能下降要小得多：8 系列有 6 个分区，而 200 系列有 8 个分区。转置内核的最终版本绝不是可以实现的最高优化水平。分块大小、每个线程的元素数量和指令优化可以提高转置和复制内核的性能。但在这项研究中，我们仅仅关注了影响最大的问题。

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 = &copy; 
 kernelName = "simple copy "; break;
 case 1: 
 kernel = &copySharedMem; 
 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 == &copy || kernel == &copySharedMem) { 
 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; 
}

可执行完整主机代码

#include <stdio.h> 
 
// 核函数将 TILE_DIM x TILE_DIM 大小的 tile 进行转置或复制
// 使用 TILE_DIM x BLOCK_ROWS 的线程块，每个线程转置 TILE_DIM/BLOCK_ROWS 个元素
// TILE_DIM 必须是 BLOCK_ROWS 的整数倍
 
#define TILE_DIM 32        // 定义 tile 的维度
#define BLOCK_ROWS 8       // 定义块中的行数

// 定义用于计时的重复次数
#define NUM_REPS 100
 __global__ void copy(float *odata, float* idata, int width, int height, int nreps) 
{ 
    // 计算当前线程在x轴上的索引
    int xIndex = blockIdx.x * TILE_DIM + threadIdx.x; 
    
    // 计算当前线程在y轴上的索引
    int yIndex = blockIdx.y * TILE_DIM + threadIdx.y; 
    
    // 计算当前线程处理的数据在一维数组中的索引位置
    int index = xIndex + width * yIndex; 
    
    // 外层循环，重复nreps次
    for (int r = 0; r < nreps; r++) { 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行处理
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将idata数组中的数据拷贝到odata数组中
            odata[index + i * width] = idata[index + i * width]; 
        } 
    } 
}
__global__ void transposeNaive(float *odata, float* idata, int width, int height, int nreps) 
{ 
    // 计算当前线程在x轴上的索引
    int xIndex = blockIdx.x * TILE_DIM + threadIdx.x; 
    
    // 计算当前线程在y轴上的索引
    int yIndex = blockIdx.y * TILE_DIM + threadIdx.y; 
    
    // 计算输入数据数组中当前线程处理的数据的位置索引
    int index_in = xIndex + width * yIndex; 
    
    // 计算输出数据数组中当前线程处理的数据的位置索引
    int index_out = yIndex + height * xIndex; 
    
    // 外层循环，重复nreps次
    for (int r = 0; r < nreps; r++) { 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行处理
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将输入数组idata中的数据拷贝到输出数组odata中，进行转置
            odata[index_out + i] = idata[index_in + i * width];
        } 
    } 
}
__global__ void transposeCoalesced(float *odata, float *idata, int width, int height, int nreps) 
{ 
    // 在共享内存中定义一个TILE_DIM x TILE_DIM的二维数组tile，用于存储矩阵块
    __shared__ float tile[TILE_DIM][TILE_DIM]; 
    
    // 计算输入数据在x轴和y轴的索引
    int xIndex = blockIdx.x * TILE_DIM + threadIdx.x; 
    int yIndex = blockIdx.y * TILE_DIM + threadIdx.y; 
    
    // 计算输入数据数组中当前线程处理的数据的位置索引
    int index_in = xIndex + (yIndex) * width; 
    
    // 重新计算xIndex和yIndex，以适应输出数据的转置
    xIndex = blockIdx.y * TILE_DIM + threadIdx.x; 
    yIndex = blockIdx.x * TILE_DIM + threadIdx.y; 
    
    // 计算输出数据数组中当前线程处理的数据的位置索引
    int index_out = xIndex + (yIndex) * height; 
    
    // 外层循环，重复nreps次
    for (int r = 0; r < nreps; r++) { 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行处理
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将输入数组idata中的数据加载到共享内存tile中
            tile[threadIdx.y + i][threadIdx.x] = idata[index_in + i * width]; 
        } 
        
        // 等待所有线程完成共享内存的数据加载
        __syncthreads(); 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行转置和写回
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将共享内存tile中的数据写入输出数组odata
            odata[index_out + i * height] = tile[threadIdx.x][threadIdx.y + i]; 
        } 
    } 
}
__global__ void transposeNoBankConflicts(float *odata, float *idata, int width, int height, int nreps) 
{ 
    // 在共享内存中定义一个TILE_DIM x TILE_DIM的二维数组tile，用于存储矩阵块
    __shared__ float tile[TILE_DIM][TILE_DIM+1];
    
    // 计算输入数据在x轴和y轴的索引
    int xIndex = blockIdx.x * TILE_DIM + threadIdx.x; 
    int yIndex = blockIdx.y * TILE_DIM + threadIdx.y; 
    
    // 计算输入数据数组中当前线程处理的数据的位置索引
    int index_in = xIndex + (yIndex) * width; 
    
    // 重新计算xIndex和yIndex，以适应输出数据的转置
    xIndex = blockIdx.y * TILE_DIM + threadIdx.x; 
    yIndex = blockIdx.x * TILE_DIM + threadIdx.y; 
    
    // 计算输出数据数组中当前线程处理的数据的位置索引
    int index_out = xIndex + (yIndex) * height; 
    
    // 外层循环，重复nreps次
    for (int r = 0; r < nreps; r++) { 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行处理
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将输入数组idata中的数据加载到共享内存tile中
            tile[threadIdx.y + i][threadIdx.x] = idata[index_in + i * width]; 
        } 
        
        // 等待所有线程完成共享内存的数据加载
        __syncthreads(); 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行转置和写回
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将共享内存tile中的数据写入输出数组odata
            odata[index_out + i * height] = tile[threadIdx.x][threadIdx.y + i]; 
        } 
    } 
}
__global__ void copySharedMem(float *odata, float *idata, int width, int height, int nreps) {
    // 在共享内存中定义一个TILE_DIM x TILE_DIM的二维数组tile
    __shared__ float tile[TILE_DIM][TILE_DIM]; 
    
    // 计算当前线程在x轴上的全局索引
    int xIndex = blockIdx.x * TILE_DIM + threadIdx.x; 
    
    // 计算当前线程在y轴上的全局索引
    int yIndex = blockIdx.y * TILE_DIM + threadIdx.y; 
    
    // 计算输入数据数组中当前线程处理的数据的位置索引
    int index = xIndex + width * yIndex; 
    
    // 外层循环，重复nreps次
    for (int r = 0; r < nreps; r++) { 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行处理
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将输入数组idata中的数据加载到共享内存tile中
            tile[threadIdx.y + i][threadIdx.x] = idata[index + i * width]; 
        } 
        
        // 等待所有线程完成共享内存的数据加载
        __syncthreads(); 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行处理
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将共享内存tile中的数据拷贝回输出数组odata中
            odata[index + i * width] = tile[threadIdx.y + i][threadIdx.x]; 
        } 
    } 
}
// 精细粒度矩阵转置核函数
__global__ void transposeFineGrained(float *odata, float *idata, int width, int height, int nreps) 
{
    // 声明一个共享内存块，用于存储 TILE_DIM x TILE_DIM 的块。
    // TILE_DIM+1 是为了避免银行冲突。
    __shared__ float block[TILE_DIM][TILE_DIM+1];

    // 计算当前线程在全局矩阵中的 x 和 y 坐标
    int xIndex = blockIdx.x * TILE_DIM + threadIdx.x;
    int yIndex = blockIdx.y * TILE_DIM + threadIdx.y;

    // 计算当前线程在输入矩阵中的线性索引
    int index = xIndex + (yIndex) * width;

    // 执行 nreps 次矩阵转置操作
    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];
        }

        // 等待所有线程都完成共享内存的读操作
        __syncthreads();
    }
}

// 粗粒度矩阵转置核函数
__global__ void transposeCoarseGrained(float *odata, float *idata, int width, int height, int nreps) 
{
    // 声明一个共享内存块，用于存储 TILE_DIM x TILE_DIM 的块。
    // TILE_DIM+1 是为了避免银行冲突。
    __shared__ float block[TILE_DIM][TILE_DIM+1];

    // 计算当前线程在输入矩阵中的 x 和 y 坐标
    int xIndex = blockIdx.x * TILE_DIM + threadIdx.x;
    int yIndex = blockIdx.y * TILE_DIM + threadIdx.y;

    // 计算当前线程在输入矩阵中的线性索引
    int index_in = xIndex + (yIndex) * width;

    // 交换 x 和 y 坐标以计算转置后的索引
    xIndex = blockIdx.y * TILE_DIM + threadIdx.x;
    yIndex = blockIdx.x * TILE_DIM + threadIdx.y;

    // 计算当前线程在输出矩阵中的线性索引
    int index_out = xIndex + (yIndex) * height;

    // 执行 nreps 次矩阵转置操作
    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];
        }

        // 等待所有线程都完成共享内存的读操作
        __syncthreads();
    }
}

__global__ void transposeDiagonal(float *odata, float *idata, int width, int height, int nreps) 
{ 
    // 在共享内存中定义一个TILE_DIM x (TILE_DIM+1)的二维数组tile
    __shared__ float tile[TILE_DIM][TILE_DIM+1]; 
    
    int blockIdx_x, blockIdx_y; 
    
    // 对角线重排序逻辑
    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; 
    } 
    
    // 计算当前线程在x轴上的全局索引
    int xIndex = blockIdx_x * TILE_DIM + threadIdx.x; 
    
    // 计算当前线程在y轴上的全局索引
    int yIndex = blockIdx_y * TILE_DIM + threadIdx.y; 
    
    // 计算输入数据数组中当前线程处理的数据的位置索引
    int index_in = xIndex + (yIndex) * width; 
    
    // 重新计算xIndex和yIndex，以适应输出数据的转置
    xIndex = blockIdx_y * TILE_DIM + threadIdx.x; 
    yIndex = blockIdx_x * TILE_DIM + threadIdx.y; 
    
    // 计算输出数据数组中当前线程处理的数据的位置索引
    int index_out = xIndex + (yIndex) * height; 
    
    // 外层循环，重复nreps次
    for (int r = 0; r < nreps; r++) { 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行处理
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将输入数组idata中的数据加载到共享内存tile中
            tile[threadIdx.y + i][threadIdx.x] = idata[index_in + i * width]; 
        } 
        
        // 等待所有线程完成共享内存的数据加载
        __syncthreads(); 
        
        // 内层循环，通过步长BLOCK_ROWS对数据进行处理
        for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) { 
            
            // 将共享内存tile中的数据写回到输出数组odata中
            odata[index_out + i * height] = tile[threadIdx.x][threadIdx.y + i]; 
        } 
    } 
}
int comparef(float* reference, float* data, int size) {
    for (int i = 0; i < size; ++i) {
        if (reference[i] != data[i]) {
            return 0;  // 不相等，返回失败
        }
    }
    return 1;  // 相等，返回成功
}
void computeTransposeGold(float* transposeGold, float* h_idata, int size_x, int size_y) {
    for (int y = 0; y < size_y; ++y) {
        for (int x = 0; x < size_x; ++x) {
            transposeGold[x * size_y + y] = h_idata[y * size_x + x];
        }
    }
}

int main(int argc, char** argv) 
{ 
    // 设置矩阵大小
    const int size_x = 2048, size_y = 2048; 
 
    // 核函数指针和描述符
    void (*kernel)(float *, float *, int, int, int); 
    const char *kernelName; 
 
    // 执行配置参数
    dim3 grid(size_x / TILE_DIM, size_y / TILE_DIM); 
    dim3 threads(TILE_DIM, BLOCK_ROWS); 
 
    // CUDA 事件，用于计时
    cudaEvent_t start, stop; 
 
    // 存储矩阵所需的内存大小
    const int mem_size = sizeof(float) * size_x * size_y; 
 
    // 分配主机内存
    float *h_idata = (float*) malloc(mem_size); 
    float *h_odata = (float*) malloc(mem_size); 
    float *transposeGold = (float *) malloc(mem_size); 
    float *gold; 
 
    // 分配设备内存
    float *d_idata, *d_odata; 
    cudaMalloc((void**) &d_idata, mem_size); 
    cudaMalloc((void**) &d_odata, mem_size); 
 
    // 初始化主机数据
    for(int i = 0; i < (size_x * size_y); ++i) 
        h_idata[i] = (float) i; 
 
    // 将主机数据复制到设备
    cudaMemcpy(d_idata, h_idata, mem_size, cudaMemcpyHostToDevice); 
 
    // 计算参考转置解
    computeTransposeGold(transposeGold, h_idata, size_x, size_y); 
 
    // 打印所有核函数的通用数据
    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"); 
 
    // 
    // 循环执行不同的核函数
    // 
    for (int k = 0; k < 8; k++) { 
        // 设置核函数指针
        switch (k) { 
            case 0: 
                kernel = &copy; 
                kernelName = "simple copy "; 
                break; 
            case 1: 
                kernel = &copySharedMem; 
                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; 
        } 
 
        // 设置参考解
        // 注意：粗粒度和细粒度的核函数不是完整的转置，因此跳过检查
        if (kernel == &copy || kernel == &copySharedMem) { 
            gold = h_idata; 
        } else if (kernel == &transposeCoarseGrained || kernel == &transposeFineGrained) { 
            gold = h_odata; 
        } else { 
            gold = transposeGold; 
        } 
 
        // 初始化事件，执行配置参数
        cudaEventCreate(&start); 
        cudaEventCreate(&stop); 
 
        // 预热以避免计时启动延迟
        kernel<<<grid, threads>>>(d_odata, d_idata, size_x, size_y, 1); 
 
        // 测量循环调用核函数的时间
        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); 
 
        // 测量核函数内部循环的时间
        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); 
 
        // 报告有效带宽
        float outerBandwidth = 
            2.0 * 1000 * mem_size / (1024 * 1024 * 1024) / (outerTime / NUM_REPS); 
        float innerBandwidth = 
            2.0 * 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); 
    } 
 
    // 清理内存
    free(h_idata); 
    free(h_odata); 
    free(transposeGold); 
    cudaFree(d_idata); 
    cudaFree(d_odata); 
    cudaEventDestroy(start); 
    cudaEventDestroy(stop); 
 
    return 0; 
}