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

前言：计算坐标过程有详图讲解哦~ヾ(◍°∇°◍)ﾉﾞ

1. 矩阵转置运算

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

$$\begin{bmatrix} 1 &2 \\ 3 &4 \\ 5 &6 \end{bmatrix}$$ A^{T}= $$\begin{bmatrix} 1 &3 &5\\ 2 &4 &6\\ \end{bmatrix}$$ A=⎣ ⎡​135​246​⎦ ⎤​AT=[12​34​56​]

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];
} 
1
2
3
4
5
6
7
8

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

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

1. row和col所用经典计算公式不再赘述，分别表示该数据所在行序号和所在列序号。
2. 从原矩阵matx的视角来看坐标，其二维展开的一维坐标为row*n+col。
3. 那么其转置矩阵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];
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();
    
    }
    
} 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

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

上述计算示意图如下所示：