前言
随着高性能计算飞速发展，异构计算已成为提升性能和效率的重要途径。特别是在图形处理单元（GPU）逐渐成为通用计算平台的背景下，GPU异构的研究和应用愈发引人注目。异构GPU计算结合了多种类型的计算资源，通过并行处理和协同工作，能够显著提高大规模数据处理、深度学习、科学计算等领域的计算能力。
GPU异构的主要平台为Cuda，基础任务是通过各种各样的软硬算法，实现矩阵乘法的优化。因为无论是模拟仿真的数值计算，还是人工智能算法的卷积、池化，归根到底还是对矩阵的处理。
而矩阵乘法的优化涉及到各种各样的矩阵拆解、计算、结合，本身对算法水平以及代码理解能力要求较高，且矩阵的拆解与计算比较抽象，过程难以模拟，此时一个帮助分析代码的“助手”就显得犹然可贵。本文主要利用通义灵码，对共享内存的矩阵乘法优化算法进行分析与拆解。
实验设备
CPU设备：i5-10210U（4核，8线程，1.60 GHz）
IDE和编译器：Visual Studio 2022/MSVC14
CUDA：12.3
优化选项：o2
通义灵码：通过Visual Studio 中安装。
具体教程如下：link
实验代码
#include <cuda_runtime.h>
#include <iostream>
#include <cstdlib>
// 定义FLOAT4宏来加载或存储4个浮点数
#define FLOAT4(x) (*((float4*)(&(x))))
// 定义OFFSET宏来计算全局内存地址
#define OFFSET(row, col, stride) ((row) * (stride) + (col))
// CUDA kernel function
__global__ void mySgemmV1Aligned(
float* __restrict__ a, float* __restrict__ b, float* __restrict__ c,
const int M, const int N, const int K)
{
const int BM = 128;
const int BN = 128;
const int BK = 8;
const int TM = 8;
const int TN = 8;
const int bx = blockIdx.x;
const int by = blockIdx.y;
const int tx = threadIdx.x;
const int ty = threadIdx.y;
const int tid = ty * blockDim.x + tx;
__shared__ float s_a[BM][BK];
__shared__ float s_b[BK][BN];
float r_c[TM][TN] = { 0.0 };
int load_a_smem_m = tid >> 1;
int load_a_smem_k = (tid & 1) << 2;
int load_b_smem_k = tid >> 5;
int load_b_smem_n = (tid & 31) << 2;
int load_a_gmem_m = by * BM + load_a_smem_m;
int load_b_gmem_n = bx * BN + load_b_smem_n;
for (int bk = 0; bk < (K + BK - 1) / BK; bk++) {
int load_a_gmem_k = bk * BK + load_a_smem_k;
int load_a_gmem_addr = OFFSET(load_a_gmem_m, load_a_gmem_k, K);
FLOAT4(s_a[load_a_smem_m][load_a_smem_k]) = FLOAT4(a[load_a_gmem_addr]);
int load_b_gmem_k = bk * BK + load_b_smem_k;
int load_b_gmem_addr = OFFSET(load_b_gmem_k, load_b_gmem_n, N);
FLOAT4(s_b[load_b_smem_k][load_b_smem_n]) = FLOAT4(b[load_b_gmem_addr]);
__syncthreads();
#pragma unroll
for (int k = 0; k < BK; k++) {
#pragma unroll
for (int m = 0; m < TM; m++) {
#pragma unroll
for (int n = 0; n < TN; n++) {
int comp_a_smem_m = ty * TM + m;
int comp_b_smem_n = tx * TN + n;
r_c[m][n] += s_a[comp_a_smem_m][k] * s_b[k][comp_b_smem_n];
}
}
}
__syncthreads();
}
#pragma unroll
for (int i = 0; i < TM; i++) {
int store_c_gmem_m = by * BM + ty * TM + i;
#pragma unroll
for (int j = 0; j < TN; j += 4) {
int store_c_gmem_n = bx * BN + tx * TN + j;
int store_c_gmem_addr = OFFSET(store_c_gmem_m, store_c_gmem_n, N);
FLOAT4(c[store_c_gmem_addr]) = FLOAT4(r_c[i][j]);
}
}
}
int main()
{
// 设定矩阵大小
const int M = 1024;
const int N = 1024;
const int K = 1024;
// 初始化主机端内存
float* h_A = new float[M * K];
float* h_B = new float[K * N];
float* h_C = new float[M * N];
// 随机初始化矩阵 A 和 B
for (int i = 0; i < M * K; ++i) h_A[i] = static_cast<float>(rand()) / static_cast<float>(RAND_MAX);
for (int i = 0; i < K * N; ++i) h_B[i] = static_cast<float>(rand()) / static_cast<float>(RAND_MAX);
// 分配设备端内存
float* d_A;
float* d_B;
float* d_C;
cudaMalloc(&d_A, M * K * sizeof(float));
cudaMalloc(&d_B, K * N * sizeof(float));
cudaMalloc(&d_C, M * N * sizeof(float));
// 将主机端数据复制到设备端
cudaMemcpy(d_A, h_A, M * K * sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(d_B, h_B, K * N * sizeof(float), cudaMemcpyHostToDevice);
// 启动CUDA内核
const int BM = 128;
const int BN = 128;
dim3 blockSize(BM, BN);
dim3 gridSize((N + BM - 1) / BM, (M + BN - 1) / BN);
cudaEvent_t start, stop;
float elapsed_time_ms = 0;
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start, 0);
mySgemmV1Aligned << <gridSize, blockSize >> > (d_A, d_B, d_C, M, N, K);
cudaEventRecord(stop, 0);
cudaEventSynchronize(stop);
cudaEventElapsedTime(&elapsed_time_ms, start, stop);
std::cout << "Matrix multiplication took " << elapsed_time_ms << " ms." << std::endl;
// 将结果从设备复制回主机
cudaMemcpy(h_C, d_C, M * N * sizeof(float), cudaMemcpyDeviceToHost);
// 释放内存
cudaFree(d_A);
cudaFree(d_B);
cudaFree(d_C);
delete[] h_A;
delete[] h_B;
delete[] h_C;
return 0;
}
基于通义灵码的代码分析
开启功能
不同于VScode，首先，需要在Visual Studio中打开通义灵码的问答框。
路径为工具->通义灵码->通义灵码问答
程序整体分析
将程序整体拖入进行分析。根据通义灵码的帮助文档，合适的提示词也是非常重要的。在这里，我使用的提示词是：
您是一位经验丰富的Cuda算法工程师，请仔细分析以下代码并对其方法与实现方式做出详细描述。
可以看到通义灵码对代码的每一个部分做了总结，并对程序整体做了介绍，使我们大体了解了程序在做什么。
程序细节分析
在Cuda优化程序中，最抽象的部分莫过于线程调度与矩阵的分块。接着，我们针对mySgemmV1Aligned()函数使用通义灵码进行单独分析。
先提问线程索引tid划分部分，为什么a和b不对称。
int load_a_smem_m = tid >> 1;
int load_a_smem_k = (tid & 1) << 2;
int load_b_smem_k = tid >> 5;
int load_b_smem_n = (tid & 31) << 2;
int load_a_gmem_m = by * BM + load_a_smem_m;
int load_b_gmem_n = bx * BN + load_b_smem_n;
得到回答
接着是外层循环
for (int bk = 0; bk < (K + BK - 1) / BK; bk++) {
int load_a_gmem_k = bk * BK + load_a_smem_k;
int load_a_gmem_addr = OFFSET(load_a_gmem_m, load_a_gmem_k, K);
FLOAT4(s_a[load_a_smem_m][load_a_smem_k]) = FLOAT4(a[load_a_gmem_addr]);
int load_b_gmem_k = bk * BK + load_b_smem_k;
int load_b_gmem_addr = OFFSET(load_b_gmem_k, load_b_gmem_n, N);
FLOAT4(s_b[load_b_smem_k][load_b_smem_n]) = FLOAT4(b[load_b_gmem_addr]);
__syncthreads();
得到回答
结合回答分析，外层循环主要起到了矩阵块切分的作用，将切分的块放到共享内存中以减少数据从全局内存中读取的时间延时。
然后是下个部分
#pragma unroll
for (int k = 0; k < BK; k++) {
#pragma unroll
for (int m = 0; m < TM; m++) {
#pragma unroll
for (int n = 0; n < TN; n++) {
int comp_a_smem_m = ty * TM + m;
int comp_b_smem_n = tx * TN + n;
r_c[m][n] += s_a[comp_a_smem_m][k] * s_b[k][comp_b_smem_n];
}
}
}
__syncthreads();
}
#pragma unroll
for (int i = 0; i < TM; i++) {
int store_c_gmem_m = by * BM + ty * TM + i;
#pragma unroll
for (int j = 0; j < TN; j += 4) {
int store_c_gmem_n = bx * BN + tx * TN + j;
int store_c_gmem_addr = OFFSET(store_c_gmem_m, store_c_gmem_n, N);
FLOAT4(c[store_c_gmem_addr]) = FLOAT4(r_c[i][j]);
}
}
}
得到回答
结合通义灵码可知，下面两步分别做了共享内存中的矩阵乘法与结果写回内存。
利用通义灵码带入具体值分析
上述部分我们从整体与局部的角度，分别利用通义灵码对代码进行了分析。但是对于矩阵，特别是二维大矩阵来说还是有些抽象，带入具体值进行分析。
如果手动拆解循环，不仅需要计算的东西很多，而且分析起来也不连贯，我们利用通义灵码进行具体的分析。我的提示词如下
我对以下循环不理解，请带入具体的值进行分析，带入不少于五组
在具体实践中，在IDE中选中代码，在问答中可以直接对选中的代码提问。
由图可见，通义灵码带入了具体的值分析了循环的过程，简洁明了。
总结
本文基于Visual Studio 2022下的通义灵码，对共享内存的矩阵乘法优化代码进行了分析，分别从整体、局部、具体值的角度，分析了代码功能及其实现方式，并演示了如何利用通义灵码进行提问与分析。
