注：这是我的线代报告，发这里算了

前言

引言

如今，线性代数在深度学习中的应用非常广泛。但是在计算机领域，使用传统CPU进行矩阵运算是非常缓慢的，因为CPU并不擅长高度并行计算。而对于擅长高度并行化计算的GPU而言，这是非常容易的，这也是为什么现代大模型的训练离不开大量显卡（计算卡）

本次报告目的

本文使用CUDA编程实现高效的矩阵操作（包括矩阵乘法、矩阵转置）

CUDA简介^[1]

CUDA 是NVIDIA 发明的一种并行计算平台和编程模型。它通过利用
图形处理器(GPU)的处理能力，可大幅提升计算性能。

编程环境

操作系统	CUDA设备	IDE
ubuntu 18 LTS	tesla p4	Nsight
windows 11	RTX 3060 laptop	CLion

基础知识

名词	中文	解释
shared memory	共享显存	容量很小，速度很快，随机读写性能好，局部访问的缓存
global memory	全局显存	容量很大，速度较慢，随机读写性能差，全局访问的缓存

矩阵转置

数学原理

我们学习的矩阵转置如下

若矩阵有A_{ij}， 则A^T=A_{ji}

但是在实际计算过程中，矩阵体积很大，只能存储在全局显存中，而矩阵转置涉及到了大量的的随机读写，导致了性能显著下降。

要解决这个问题，就需要引入共享显存，但是共享显存的体积很小且只能局部访问，无法存储整个矩阵。所以，我们可以使用分块矩阵来优化这个过程。将子矩阵复制到共享显存，再使用共享显存完成子矩阵的转置，然后再复制回全局显存。

分块矩阵转置： \begin{bmatrix} A&B\\ C&D \end{bmatrix}^T = \begin{bmatrix} A^T&C^T\\ B^T&D^T \end{bmatrix}

CUDA代码实现

__global__ void kernel_transpose(DTYPE *input, DTYPE *output, int num_rows, int num_cols)
{
    __shared__ DTYPE buffer[BLOCK_SIZE][BLOCK_SIZE];//创建共享显存
    int column=blockIdx.x*BLOCK_SIZE+threadIdx.x;
    int row=blockIdx.y*BLOCK_SIZE+threadIdx.y;
    if(column>=num_cols||row>=num_rows){//如果行数或列数超过了矩阵大小，则关闭这个线程
        return;
    }
    buffer[threadIdx.y][threadIdx.x]=input[row*num_cols+column];//转置并存入共享显存
    __syncthreads();//等待所有线程完成该操作
    output[column*num_rows+row]=buffer[threadIdx.y][threadIdx.x];//将共享显存的数据复制回全局显存
}

矩阵乘法

数学原理

我们学习的矩阵乘法如下

\begin{aligned} 有(M)_{ij},(N)_{jm},P=M \times N\\ 则(P)_{im}=\sum_{k=1}^{j}M_{ik}N_{kj} \end{aligned}

这样的确是最简单的算法，使用Python语言的列表推导式能在一行内完成这个操作

1 2	def mul(m,n): return [[sum([m[i][k]*n[k][j] for k in range(len(n))]) for j in range(len(n[0]))] for i in range(len(m))]

但是由于与上文相同的原因，大量的随机写入导致性能大幅下降。因此，必须使用分块矩阵乘法。数学原理与上文类似，不过是将每一个数字换成了矩阵，对于每一个小矩阵，放在共享显存中完成运算，再复制回全局显存。这样大大减少了对全局显存的随机读取，大幅提高了计算效率。

CUDA代码实现

__global__ void kernel_matrix_multiply_block(float *M, float *N, float *P, int M_rows, int M_cols, int N_rows, int N_cols)
{
    __shared__ float M_shared[BLOCK_SIZE][BLOCK_SIZE];//创建M的分块矩阵共享显存
    __shared__ float N_shared[BLOCK_SIZE][BLOCK_SIZE];//创建N的分块矩阵共享显存

    int bx = blockIdx.x;
    int by = blockIdx.y;
    int tx = threadIdx.x;
    int ty = threadIdx.y;

    int row = by * BLOCK_SIZE + ty;
    int col = bx * BLOCK_SIZE + tx;

    if (row >= M_rows || col >= N_cols) {
        return;
    }

    float Pvalue = 0;
    for (int m = 0; m < (M_cols - 1) / BLOCK_SIZE+1; m++) {
        //将分块矩阵（子矩阵）复制到共享显存中
        if (row < M_rows && m * BLOCK_SIZE + tx < M_cols) {
            M_shared[ty][tx] = M[row * M_cols + m * BLOCK_SIZE + tx];
        } else {
            M_shared[ty][tx] = 0;
        }
        if (col < N_cols && m * BLOCK_SIZE + ty < N_rows) {
            N_shared[ty][tx] = N[(m * BLOCK_SIZE + ty) * N_cols + col];
        } else {
            N_shared[ty][tx] = 0;
        }
        __syncthreads();
        //在共享显存中完成矩阵乘法
        for (int k = 0; k < BLOCK_SIZE; ++k) {
            Pvalue += M_shared[ty][k] * N_shared[k][tx];
        }
        __syncthreads();
    }
    //将乘法结果写回全局显存
    P[row * N_cols + col] = Pvalue;
    
}