矩阵转置及优化——CUDA on ARM Platform夏令营 Day03 课后作业

Lotus莲

已于 2022-07-15 16:36:19 修改

阅读量614

点赞数 1

分类专栏： 2022 CUDA Summer Cap 文章标签：矩阵 arm 线性代数

于 2022-07-07 00:34:31 首次发布

本文链接：https://blog.csdn.net/weixin_48266461/article/details/125647836

版权

2022 CUDA Summer Cap 专栏收录该内容

2 篇文章 0 订阅

订阅专栏

1. 矩阵转置运算

大小为m×n的矩阵 $A$ ，其转置矩阵 $A^{T}$ 的大小为n×m，且有矩阵 $A$ 的第 $i$ 行第 $j$ 列的元素 $a (i, j)$ 与矩阵 $A^{T}$ 的第 $j$ 行第 $i$ 列的元素 $a^{T}(j,i)$ 相等。
举例如下：
$A=\begin{bmatrix} 1 &2 \\ 3 &4 \\ 5 &6 \end{bmatrix} A^{T}=\begin{bmatrix} 1 &3 &5\\ 2 &4 &6\\ \end{bmatrix}$

2. 矩阵转置GPU代码

先上代码，解释在后~

__global__ void gpu_matrix_trans(int *matx, int *matxT, int m, int n)
{ 
    int row = blockIdx.y * blockDim.y + threadIdx.y; 
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    if( row < m && col < n) 
        matxT[col*m+row] = matx[row*n+col];
}

这里是将二维矩阵展开为一维进行表示（申请的内存为连续的一维形式）。

cudaMalloc((void **) &matx, sizeof(int)*m*n)

row和col所用经典计算公式不再赘述，分别表示该数据所在行序号和所在列序号。
从原矩阵matx的视角来看坐标，其二维展开的一维坐标为row*n+col。
那么其转置矩阵matxT对应的一维坐标为col*m+row。

坐标计算示意图如下所示：
在这里插入图片描述

3. 运用共享内存进行矩阵转置优化

3.1 共享内存

相较于上述数据所使用的全局内存（Global memory），使用共享内存（Shared memory）可以对矩阵转置进行进一步优化。
共享内存位于GPU芯片上，相较位于DRAM中的全局内存拥有更低的延迟，其存取速度仅次于最快的寄存器。
不同的存储单元
（上图来自夏令营PPT）
但共享内存大小有限，数据量较大时需要对数据进行划分处理，否则会限制活动warp的数量。实现过程将在后续代码中体现。

3.2 bank conflict

同一block中的线程都可以访问共享内存，因此可能出现同时有很多线程访问共享内存中的数据。为了克服此问题，共享内存被划分为32个大小相等的内存块bank。
但同一warp中的多个线程同时访问同一个bank时就会发生冲突，降低效率。（但若访问同一地址，则会触发广播机制，不会发生冲突）
在这里插入图片描述

（上图来自夏令营PPT）
如上图所示，为避免冲突，可以在原申请内存的基础上增加一列，所以这里需要在申请共享内存时给BLOCK_SIZE+1。

__shared__ int tile_p[BLOCK_SIZE][BLOCK_SIZE+1];

3.3 实现代码

#define BLOCK_SIZE 16

__global__ void gpu_matrix_trans_shared(int *d_matx,int *d_matxTs, int m, int n)
{ 
    //申请共享内存
    __shared__ int tile_p[BLOCK_SIZE][BLOCK_SIZE+1];
    
    //转置矩阵元素的位置
    int row = blockIdx.y * blockDim.y + threadIdx.y; 
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    
    //循环“小块”转置
    int idx;
    for(int sub=0; sub<(n/BLOCK_SIZE+1); sub++){
        //将数据copy到shared memory里面
        idx = row*n + BLOCK_SIZE*sub + threadIdx.x;
        if(row < m && (BLOCK_SIZE*sub + threadIdx.x) < n)
            tile_p[threadIdx.y][threadIdx.x] = d_matx[idx];
            
        //转置结果
        __syncthreads();
        idx = (BLOCK_SIZE*sub + threadIdx.x)*m + row;
        if(row < m && (BLOCK_SIZE*sub + threadIdx.x) < n)
            d_matxTs[idx] = tile_p[threadIdx.y][threadIdx.x];
        __syncthreads();
    
    }
    
}

首先使用__shared__申请共享内存块。
同样地得到row和col。
使用for循环每一“小块”每一“小块”地处理，这里的循环次数即为“小块”在n方向滑动处理完所有数据所需的次数，所以为sub<(n/BLOCK_SIZE+1)。（C语言默认int类型相除为向下取整）
下面是很关键的坐标计算，需要找准每一小块和矩阵、矩阵和其转置矩阵间的坐标关系。对于小块tile_p[threadIdx.y][threadIdx.x]其对应的矩阵matx的展开的一维坐标为row* n + BLOCK_SIZE*sub + threadIdx.x。这里还要注意使用正确约束条件进行越界判断，保证row < m && (BLOCK_SIZE*sub + threadIdx.x) < n，否则会带来问题。
与上述原矩阵坐标对应的转置矩阵matxT对应的一维坐标为(BLOCK_SIZE* sub + threadIdx.x)* m + row，判断条件为row < m && (BLOCK_SIZE*sub + threadIdx.x。
还有记得使用__syncthreads();完成同步！

上述计算示意图如下所示：
在这里插入图片描述