CUDA学习之路[11]：矩阵乘法计算详解
写在前边
写了这么久，终于写到了矩阵运算，后续的代码普遍都比较长，比较综合，完整的代码仓库在这里，希望对大家有所帮助和收获！
什么是矩阵乘法
无论是卷积层、全连接层还是注意力计算，最终都会归因到矩阵运算。 无论是C语言还是CUDA亦或者是numpy、pytorch，都有一大堆针对矩阵运算来进行优化的库，这也是我们AI Infra真正意义上的起点，跟随着我的视角，我们来完整的体会一下是如何优化的。
我们先来看一下CPU版本的矩阵运算：
void sgemm(const float* A, const float* B, float* C, int M, int N, int K) { for (int i = 0; i < M; i++) { for (int j = 0; j < N; j++) { float sum = 0.0f; for (int k = 0; k < K; k++) { sum += A[i * K + k] * B[k * N + j]; } C[i * N + j] = sum; } }}先分析一下这个算法的时间复杂度，三个for循环，是。
好，我们已经学过了CUDA、Triton、Pytorch。他们分别是从线程、块、张量的角度来进行分析的。
那矩阵乘法应该如何优化呢？ 所有优化的核心都只有两点：
- 优化访存
- 优化计算
Pytorch
在Pytorch中，有非常多的计算矩阵乘法的方式，最常见的是以下几种：
@与*的区别
对于两个相同尺寸的张量来说，*运算和@或者matmul运算是不同的。
*运算符表示各个元素直接相乘，例如a*b，它的结果直接就是对应元素相乘。
但是@运算符表示的是矩阵相乘
a = torch.tensor([1,2])b = torch.tensor([2,3])
print(a*b)print(torch.matmul(a, b))print(a @ b)
# tensor([2, 6])# tensor(8)# tensor(8)矩阵相乘的若干写法
除了@和matmul，在Pytorch中，还有很多中矩阵相乘的写法，对于两个向量来说，dot、@、matmul这三种写法都是等价的，其中dot只接受一维向量计算。
# @ 最常见的运算符v1 = torch.tensor([1.0, 2.0, 3.0])v2 = torch.tensor([4.0, 5.0, 6.0])
# 以下三种写法完全等价print(torch.matmul(v1, v2))print(v1 @ v2)print(torch.dot(v1, v2)) # dot 只接受 1D其中对于纯二维矩阵来说，以下三种写法是比较常见的，mm表示纯二维矩阵乘。
# 纯二维矩阵乘法a = torch.randn(100, 200)b = torch.randn(200, 300)c = torch.mm(a, b)d = torch.matmul(a, b)e = a @ b对于批量三维矩阵乘法来说，bmm是更加高效的选择：
# 批量三维矩阵乘法 (Batch, seq, Feature)batch_x = torch.randn(10, 16, 64)batch_w = torch.randn(10, 64, 128)y = torch.bmm(batch_x, batch_w)还有能够显示指定的方案：
# einsum爱因斯坦求和，通过下标字符串实现复杂的乘法、转置以及缩并Q = torch.randn(2, 8, 10, 64) # (batch, heads, len_q, d_k)K = torch.randn(2, 8, 10, 64)
# 计算 QK^Tattn = torch.einsum('b h q d, b h k d -> b h q k', Q, K) # (2,8,10,10)对于mm、bmm、dot来说，这些是最底层的算子，只接受严格的对应维度的张量，它们的低层是基于Lua的Torch7。而Matmul更加广泛，支持多种维度的张量相乘，并且同样支持广播机制。例如(10,3,4)可以和(4, 5)的张量相乘，会把后者当成(1,4,5)并且张量到(10,4,5)来计算。而@则是matmul的语法糖。但是matmul在真正执行前需要跑一大段C++的推断，多了一些检查策略。并且广播机制真的很容易引发报错。但是高纬度张量，还不得不使用这个，要不然就得疯狂使用view和transpose。至于einsum，对于高维张量来了说，这个方法我感觉是最清晰的。
如果没有@，对于A，B，C这种连乘，就得写torch.matmul(torch.matmul(A, B), C)， 有了@直接写 A @ B @ C，优雅~
Triton
Triton的核心分为两个模块，Pytorch封装和Kernel代码。
先写Pytorch封装端，做好维度检查与内存申请：
def matmul(a: torch.Tensor, b: torch.Tensor, BLOCK_SIZE_M=64, BLOCK_SIZE_N=64, BLOCK_SIZE_K=32) -> torch.Tensor: assert a.ndim == 2 and b.ndim == 2, "只支持二维矩阵" assert a.shape[1] == b.shape[0], f"维度不匹配：{a.shape} 和 {b.shape}" M, K = a.shape K, N = b.shape
# 2. 申请C的内存 c = torch.empty((M, N), device=a.device, dtype=a.dtype)在做好准备工作之后，需要把大矩阵切分成为很多个Block块，每一个块都让GPU去算：
# M 方向的块数 = ceil(M / BLOCK_SIZE_M) # N 方向的块数 = ceil(N / BLOCK_SIZE_N) # 这是一个匿名函数，meta 里面会传入 BLOCK_SIZE 的实际值 grid = lambda meta: ( triton.cdiv(M, meta["BLOCK_SIZE_M"]) * triton.cdiv(N, meta["BLOCK_SIZE_N"]), )做好grid切分之后，启动Kernel，这里要注意一个细节，通常在传入参数的时候，需要以下内容：
- 张量指针
- 张量形状
- 内存步长
之前讲过，在底层，矩阵是一维数组，并且有可能是非连续，因此显式传入步长会更有效。
matmul_kernel_v0[grid]( a, b, c, # 传张量（Triton会自动提取它们的显存首地址指针） M, N, K, # 传矩阵真实大小 a.stride(0), a.stride(1), # A 的行步长、列步长 b.stride(0), b.stride(1), # B 的行步长、列步长 c.stride(0), c.stride(1), # C 的行步长、列步长 BLOCK_SIZE_M=BLOCK_SIZE_M, # 传编译期常量 BLOCK_SIZE_N=BLOCK_SIZE_N, BLOCK_SIZE_K=BLOCK_SIZE_K, ) return c在写好了启动函数之后，接下来就是针对一个具体的BLOCK来写Kernel函数。
常见的Triton思路还是：
- 计算pid
- 计算offset和mask
- load
- 计算
- store
具体来说：
# 查看每一个block的pid pid = tl.program_id(axis=0)
# 算一下 C 矩阵一行有多少个块？ num_pid_n = tl.cdiv(N, BLOCK_SIZE_N)
# 把一维块变成二维坐标 (行块号, 列块号) pid_m = pid // num_pid_n pid_n = pid % num_pid_n
# 计算当前pid具体要负责原矩阵 C 的哪些行号和列号？ # tl.arange(0, 64) 生成 [0, 1, ... 63] off_m = pid_m * BLOCK_SIZE_M + tl.arange(0, BLOCK_SIZE_M) # 比如 [64, 65... 127] off_n = pid_n * BLOCK_SIZE_N + tl.arange(0, BLOCK_SIZE_N) # 比如 [0, 1... 63]
# K 是我们要循环的维度，每次处理一个小段 off_k = tl.arange(0, BLOCK_SIZE_K)做好数据处理后，需要算出对应的显存地址，这里使用到了广播机制：
# off_m[:, None] 变成 64x1 的列向量，乘以行步长 # off_k[None, :] 变成 1x32 的行向量，乘以列步长 a_ptrs = a_ptr + off_m[:, None] * stride_am + off_k[None, :] * stride_ak b_ptrs = b_ptr + off_k[:, None] * stride_bk + off_n[None, :] * stride_bn
# 寄存器，全部填 0，用来一会做累加,为了不掉精度，这里用 float32 acc = tl.zeros((BLOCK_SIZE_M, BLOCK_SIZE_N), dtype=tl.float32)之后沿着K的维度来进行计算：
for k in range(0, tl.cdiv(K, BLOCK_SIZE_K)):
# 边界检查：如果 K 不是 32 的整数倍，最后一步可能会越界 k_remain = K - k * BLOCK_SIZE_K
# 生成Mask。如果坐标越界了，那对应的地方就是 False a_mask = (off_m[:, None] < M) & (off_k[None, :] < k_remain) b_mask = (off_k[:, None] < k_remain) & (off_n[None, :] < N)
# 从显存搬砖到SRAM。遇到越界的地方 (Mask=False)，就用 0.0 填补 a = tl.load(a_ptrs, mask=a_mask, other=0.0) b = tl.load(b_ptrs, mask=b_mask, other=0.0)
# 算这一个小块的矩阵乘法，结果加到acc 上 acc += tl.dot(a, b, allow_tf32=False)
# 算完这一段了，把 A 和 B 的指针往 K 的方向推进一步，去拿下一段数据 a_ptrs += BLOCK_SIZE_K * stride_ak b_ptrs += BLOCK_SIZE_K * stride_bk最终写回到C中：
c_ptrs = c_ptr + off_m[:, None] * stride_cm + off_n[None, :] * stride_cn
# 15. C 的边界保护遮罩 c_mask = (off_m[:, None] < M) & (off_n[None, :] < N)
# 把 float32 的寄存器结果，转回和矩阵 C 一样的数据类型 (比如 float16) tl.store(c_ptrs, acc.to(c_ptr.dtype.element_ty), mask=c_mask)至于V1版本，实际上可以进行Group优化，那么在这里，大家可以参考Triton官方的版本，写的也是非常不错的，我就不画蛇添足了。
CUDA
V0 Naive
在CPU的实现中，两个外层循环负责遍历输出矩阵C的所有元素，最内层循环负责进行点积运算。
计算C[i][j]的过程是顺序的，先算C[0][0], 再算C[0][1]…
GPU的思路是什么呢？把个输出元素的计算任务，一次性分配给个线程。每个线程负责计算对应的C[row][col]，各自独立地从A取一整行，从B取一整列，完成K次乘加。
矩阵乘法中输出元素之间是没有数据依赖的，这就是典型的element wise并行模式。GPU Naive正是利用了这一点，用并行机制消除了两层外层循环。
最基本的V0代码如下：
__global__ void sgemm_v0(const float* A, const float* B, float* C, int M, int N, int K) { int row = blockIdx.y * blockDim.y + threadIdx.y; int col = blockIdx.x * blockDim.x + threadIdx.x; if (row < M && col < N) { float sum = 0.f; for (int k = 0; k < K; ++k) sum += A[row * K + k] * B[k * N + col]; C[row * N + col] = sum; }}V1 Tile版本
V0版本解决了两层for循环的问题，那么它如何进行性能提升呢？
大家还记得GPU上分为了HBM、Shared Memory、Register至少三级存储吗？ 在V0这种模式下，同一个数据会被不同的线程反复读取，这种读取会占用大量的带宽，并且数据也没有被进行复用。
举个例子：
计算C[0][0]线程读取A的第0行，读取B的第0列
计算C[0][1]线程读取A的第0行，读取B的第1列理论上，A的一行会被N个线程各读一遍，B的一列也会被M个线程各读一遍，这导致了巨大的冗余访存，所以这也是我们优化的方向。
一个很自然的想法是，既然A的行和B的列会被多次用到，如果我们能够先搬一部分到更快的共享内存上，各自的线程先计算一部分，这样就能够大幅度减少全局内存访问。
请看如下演示：
- 把输出矩阵分成若干[BM，BN]的块
- 计算过程中，A和B沿着K轴分成BK大小的块
- 每次迭代把A[BM][BK]和B[BK][BN]加载到共享内存中
- 线程块内所有线程从共享内存中读数据，完成局部乘加
- 循环到下一个K分块，直到遍历完整个K
每个A、B中的元素被加载到共享内存后，可以由线程块内个线程共用，全局访存量下降为原来的(1/(BM * BN))量级（粗略估计），缓解带宽瓶颈。
具体的Kernel代码如下：
__global__ void sgemm_v1(const float* A, const float* B, float* C, int M, int N, int K) { constexpr int BM = 32, BN = 32, BK = 32; __shared__ float SA[BM][BK], SB[BK][BN];
// 单个Block内部的偏移量 int ty = threadIdx.y; int tx = threadIdx.x;
// C[row][col] int row = blockIdx.y * BM + ty; int col = blockIdx.x * BN + tx;
float sum = 0.0f;
for (int tile = 0; tile < K; tile += BK) { int ga = tile + tx; SA[ty][tx] = (row < M && ga < K) ? A[row * K + ga] : 0.0f; int gb = tile + ty; SB[ty][tx] = (gb < K && col < N) ? B[gb * N + col] : 0.0f; __syncthreads();
#pragma unroll for (int k = 0; k < BK; k++) { sum += SA[ty][k] * SB[k][tx]; } __syncthreads(); }
if (row < M && col < N) { C[row * N + col] = sum; }}性能对比如下：
v0-naive 0.4501 ms 4.77 TFLOPS refv1-smem 0.3552 ms 6.05 TFLOPS 0.00e+00V1 1D Thread
在V1中，每个线程负责一个32*32的输出子块，线程块也是32*32，恰好一个线程负责一个输出元素。
在一次Tile迭代中，只做了2次global load，进行了32次shared memory load，还有32次FMA，然后就开始sync。
线程的大部分时间都在等待同步、等待load完成，计算密度太低。
优化访存我们已经在v1做过了，那是不是接下来就要尝试着优化计算了呢？
实际上，1D Tiling的核心思想是让一个线程负责多个元素，而不是一个。 原本1024个线程对应1024个位置，但是实际上我们可以让一个线程负责多个位置的计算。
线程(ty=0,tx=0)可以负责row[0, 15]， 线程(ty=1,tx=0)负责row[16, 31]…
所有的(ty, tx)线程的列 col = blockCol + tx。
这就是我们这个方法的核心，通过增加计算密度，还overlap访存效率的低下。
示意图如下：
具体代码如下：
template <int BM = 128, int BN = 64, int BK = 4, int TM = 16>__global__ void sgemm_v1_1(const float* __restrict__ A, const float* __restrict__ B, float* __restrict__ C, int M, int N, int K) {
constexpr int NT = (BM / TM) * BN; __shared__ float SA[BM][BK]; __shared__ float SB[BK][BN];
int tx = threadIdx.x; int ty = threadIdx.y;
int col = blockIdx.x * BN + tx;
float sum[TM] = {0.0f};
int tid = ty * blockDim.x + tx;
for (int tile = 0; tile < K; tile += BK) {
#pragma unroll for (int i = tid; i < BM * BK; i += NT) { int load_a_row = i / BK; int load_a_col = i % BK; int ga_row = blockIdx.y * BM + load_a_row; int ga_col = tile + load_a_col; SA[load_a_row][load_a_col] = (ga_row < M && ga_col < K) ? A[ga_row * K + ga_col] : 0.0f; }
#pragma unroll for (int i = tid; i < BK * BN; i += NT) { int load_b_row = i / BN; int load_b_col = i % BN; int gb_row = tile + load_b_row; int gb_col = blockIdx.x * BN + load_b_col; SB[load_b_row][load_b_col] = (gb_row < K && gb_col < N) ? B[gb_row * N + gb_col] : 0.0f; }
__syncthreads();
#pragma unroll for (int k = 0; k < BK; k++) { float b_val = SB[k][tx];
#pragma unroll for (int m = 0; m < TM; m++) { sum[m] += SA[ty * TM + m][k] * b_val; } } __syncthreads(); }
#pragma unroll for (int m = 0; m < TM; m++) { int row = blockIdx.y * BM + ty * TM + m;
if (row < M && col < N) { C[row * N + col] = sum[m]; } }}核心加载数据其实也没有发生明显变化，最为核心的还是在单个线程的计算上，之前是一行一列。
而现在一列要对应16行，例如选出一个B的元素，把它放在寄存器中，之后就和A的16行来进行计算，这样就进一步提升了线程的计算量。
V1 2D Thread
既然我们已经把使用单个线程把B的元素复用了，那实际上A的元素也可以明显复用，这样对于一个线程，复用率会更多，有点难以想象？看一下这个动画：
在1D的时候，是两个向量进行乘加，最终是一个标量，也就是内积：
但是矩阵乘法也可以等价地写成外积的形式。
在此感谢大佬的仔细讲解
外积把A变成列的形式，把B变成行的形式，这样同时取A的列，B的行，就完成了一次外积，然后把所有的A列和所有的B行算出来的矩阵加起来，这就是外积的形式。
例如：
A = [[1,1], [3,4]]B = [[2,2], [3,4]]# 内积形式C[0][1] = 1*2 + 1*3 = 5
# 外积形式C[0] = [1,3]^T @ [2,2] = [[2,2], [6,6]]根据这个外积，我们就可以完整的把2D Tiling完整的写出来：
template <int BM = 128, int BN = 128, int BK = 16, int TM = 8, int TN = 8>__global__ void sgemm_v1_2(const float* __restrict__ A, const float* __restrict__ B, float* __restrict__ C, int M, int N, int K) { constexpr int NT = (BM / TM) * (BN / TN); __shared__ float SA[BM][BK]; __shared__ float SB[BK][BN];
int tx = threadIdx.x; int ty = threadIdx.y; int tid = ty * blockDim.x + tx;
float sum[TM][TN] = {{0.f}};
for (int tile = 0; tile < (K + BK - 1) / BK; ++tile) { int gk = tile * BK;
#pragma unroll for (int i = tid; i < BM * BK; i += NT) { int load_a_row = i / BK; int load_a_col = i % BK; int ga_row = blockIdx.y * BM + load_a_row; int ga_col = gk + load_a_col; SA[load_a_row][load_a_col] = (ga_row < M && ga_col < K) ? A[ga_row * K + ga_col] : 0.f; }
#pragma unroll for (int i = tid; i < BK * BN; i += NT) { int load_b_row = i / BN; int load_b_col = i % BN; int gb_row = gk + load_b_row; int gb_col = blockIdx.x * BN + load_b_col; SB[load_b_row][load_b_col] = (gb_row < K && gb_col < N) ? B[gb_row * N + gb_col] : 0.f; } __syncthreads();
#pragma unroll for (int k = 0; k < BK; ++k) { float a_reg[TM], b_reg[TN]; #pragma unroll for (int m = 0; m < TM; ++m) a_reg[m] = SA[ty * TM + m][k]; #pragma unroll for (int n = 0; n < TN; ++n) b_reg[n] = SB[k][tx * TN + n]; #pragma unroll for (int m = 0; m < TM; ++m) #pragma unroll for (int n = 0; n < TN; ++n) sum[m][n] += a_reg[m] * b_reg[n]; } __syncthreads(); }
#pragma unroll for (int m = 0; m < TM; ++m) { int row = blockIdx.y * BM + ty * TM + m; if (row >= M) continue; #pragma unroll for (int n = 0; n < TN; ++n) { int col = blockIdx.x * BN + tx * TN + n; if (col < N) C[row * N + col] = sum[m][n]; } }}到此为止的性能对比情况：
kernel time TFLOPS vs ref v0-naive 0.4817 ms 4.46 TFLOPS ref v1-smem 0.3805 ms 5.64 TFLOPS 0.00e+00 v1-1Dtiling 0.2414 ms 8.90 TFLOPS 0.00e+00 v1-2Dtiling 0.1314 ms 16.34 TFLOPS 0.00e+00V2 Tile
关于计算，我们已经相比较原来提升了约十几倍，那么，还有什么地方可以提升呢？
之前是我们把数据从HMB搬到了Shared Memory，对于“搬运”这个动作来说，实际上也可以提升，例如——————————————
大声告诉我是什么！
对，就是向量化，把所有和数据搬运相关的，都改成向量化即可：
template <int BM = 64, int BN = 64, int BK = 32, int TM = 8, int TN = 4>__global__ void sgemm_v2(const float* __restrict__ A, const float* __restrict__ B, float* __restrict__ C, int M, int N, int K) {
static_assert(BK % 4 == 0, "BK must be multiple of 4 for float4"); static_assert(BN % 4 == 0, "BN must be multiple of 4 for float4"); static_assert(TN % 4 == 0, "TN must be multiple of 4 for float4");
constexpr int NT = (BM / TM) * (BN / TN); __shared__ float SA[BM][BK]; __shared__ float SB[BK][BN];
int tx = threadIdx.x; int ty = threadIdx.y; int tid = ty * blockDim.x + tx;
float sum[TM][TN] = {{0.f}};
for (int tile = 0; tile < (K + BK - 1) / BK; ++tile) { int gk = tile * BK;
#pragma unroll for (int i = tid; i < BM * BK / 4; i += NT) { int load_a_row = i / (BK / 4); int load_a_col = (i % (BK / 4)) * 4; int ga_row = blockIdx.y * BM + load_a_row; int ga_col = gk + load_a_col; if (ga_row < M && ga_col < K) *(float4*)&SA[load_a_row][load_a_col] = *(const float4*)&A[ga_row * K + ga_col]; else *(float4*)&SA[load_a_row][load_a_col] = make_float4(0.f, 0.f, 0.f, 0.f); }
#pragma unroll for (int i = tid; i < BK * BN / 4; i += NT) { int load_b_row = i / (BN / 4); int load_b_col = (i % (BN / 4)) * 4; int gb_row = gk + load_b_row; int gb_col = blockIdx.x * BN + load_b_col; if (gb_row < K && gb_col < N) *(float4*)&SB[load_b_row][load_b_col] = *(const float4*)&B[gb_row * N + gb_col]; else *(float4*)&SB[load_b_row][load_b_col] = make_float4(0.f, 0.f, 0.f, 0.f); } __syncthreads();
#pragma unroll for (int k = 0; k < BK; ++k) { #pragma unroll for (int n = 0; n < TN; n += 4) { float4 b_val = *(float4*)&SB[k][tx * TN + n]; #pragma unroll for (int m = 0; m < TM; ++m) { float a_val = SA[ty * TM + m][k]; sum[m][n + 0] += a_val * b_val.x; sum[m][n + 1] += a_val * b_val.y; sum[m][n + 2] += a_val * b_val.z; sum[m][n + 3] += a_val * b_val.w; } } } __syncthreads(); }
#pragma unroll for (int m = 0; m < TM; ++m) { int row = blockIdx.y * BM + ty * TM + m; if (row >= M) continue; #pragma unroll for (int n = 0; n < TN; n += 4) { int col = blockIdx.x * BN + tx * TN + n; if (col < N) *(float4*)&C[row * N + col] = make_float4(sum[m][n], sum[m][n + 1], sum[m][n + 2], sum[m][n + 3]); } }}V3 双缓冲
无论是之前的Norm还是其他系列，我都建议，只要是稍微复杂一些的load和store，涉及到向量化的，都可以把他封装起来，这样效果会比较好。
已经有了向量化和2D Tile，接下来还有什么优化空间呢？
我先卖个关子，可以先看一下最终优化的效果：
CUDA SGEMM — v0 to v6 optimization journey GPU: NVIDIA GeForce RTX 4090 D | SMs: 114 | sm_89--- SGEMM M=N=K=4096 (2.1 GFLOPs) --- kernel time TFLOPS vs ref v0-naive 0.4817 ms 4.46 TFLOPS ref v1-smem 0.3804 ms 5.65 TFLOPS 0.00e+00 v1-1Dtiling 0.1672 ms 12.84 TFLOPS 0.00e+00 v1-2Dtiling 0.1148 ms 18.70 TFLOPS 0.00e+00 v2-tile 0.0735 ms 29.23 TFLOPS 0.00e+00 v3-dbuf 0.0687 ms 31.25 TFLOPS 0.00e+00 cuBLAS(FP32) 0.0650 ms 33.05 TFLOPS 1.57e-03最终的V3版本大概只和cuBLAS有0.002ms的性能差距，甚至cuBLAS还偷偷用到了其他硬件模块进行加速。
那这个效果可以说是非常不错。
是怎么实现的呢？ 在V2中，sA布局是行优先，和Global A的布局是一致的，但是在访问的时候：
float a_val = SA[ty * TM + m][k]// 线程ty=0，m=0读SA[0][k]// 线程ty=0, m=1读SA[1][k]我们发现，SA[0][K]和SA[1][k]在内存中相距BK个float，这是不连续的，因此没办法做到float4读取，只能在B中float，而在A中无法这样做。
怎么解决呢？——转置存储 只需要把A的一行散射到SA的4个不同的行即可，这样就访问的是连续的行号：
*(float4*)(av+m) = *(const float4*)&SA[rb][k][ty * TM + m];除了转置之外，Double Buffering也是一个非常关键的思路。
什么意思，我们具体来看：
V2的执行流：tile 0: load SA, SB -> sync -> compute -> synctile 1: load SA, SB -> sync -> compute -> sync...如果学过体系结构的同学已经看出来可以优化的点了：Pipeline。
实际上Load和Compute是严格拆串行的，但是Tile0和Tile1之间的关联又没有那么明显，因此我们可以尝试pipeline：
V3的执行流：tile i : global -> register 加载tile i-1: shared -> register 计算...这两个同时发生，因此访存被进一步隐藏。
最终的性能想必大家也看到了，基本在cuBLAS的95%左右，当然，我还需要微调模板参数，来让他性能达到最优，这里我就先按下不表。
到此为止，和矩阵乘法的CUDA相关的内容就写的差不多了，代码仓库在文章开头，大家可以自行下载，多多Star，感谢大家的支持！ 希望对大家有所帮助。
支持与分享
如果这篇文章对你有帮助，欢迎分享给更多人或赞助支持！
音乐
暂未播放