写一个 RTX 5080 上的 cuda gemm fp16

1. cpu 计算 fp16 四则运算

由于会用到cpu 的gemm 与 gpu gemm 的对比验证，所以，这里稍微解释一下 cpu 计算fp16 gemm 的过程。这里为了简化理解，cpu 中不使用 avx 相关的 fp16 运算器，而是直接使用 cpu 原先的 ALU 功能。这里使用一个示例来做这件事情。

1.1. 源码编译运行

hello_fp16.cu


#include <stdio.h>
#include "cuda_fp16.h"int main()
{half x = half(3.333);half y = half(7.777);half z = half(0.0);z = x*y;printf("sizeof(half) = %ld x = %f \n", sizeof(x), float(z));return 0;
}

编译运行：

nvcc -g --gpu-architecture=sm_120 hello_fp16.cu -o hello_fp16

1.2. 调试追踪 fp16 的相关功能

这里有两个目标：

一个是类型转换，怎么样得到一个 fp16 的变量值；

一个是 fp16 类型变量之间的乘法（四则运算）。

现在看一下其中的 half(3.333) 的执行，通过使用 step，经历如下几个断点：

在 408 行时，使用gdb ）s 会跳到下图代码 549 行处：

继续使用 (cuda-gdb) s 会跳到下图：

然后使用 (cuda-gdb) next ，会经历上图代码的主体逻辑，也就是一些位运算的那些逻辑。

结论便是，cuda 程序对 cpu 的 half(3.333) 使用了 cpu 软件算法模拟了这个转换过程，将 double 类型转换成 fp16，即 half 类型 。

同时接下来会发现，两个 half 类型的变量做乘法运算，会先将两个 half 转成 float，也是通过类似的软件模拟的转换方式，然后使用 cpu 的 float 乘法指令计算乘积，最后将 float 类型的乘积再转回 half 类型，存入 half 类型的变量内存中。

接下来调试 half 的乘法运算符 * ：

在执行到 z = x*y; 时，使用(cuda-gdb) step,会跳进half 类型的乘法运算符 * 的实现代码中，这里使用了 cpp 的重载功能(Operator Overloading) ，对运算符 * 做了重新实现：

可以看到，operator * 重载时，函数体中调用了 __hmul(...) 来实现具体功能。

接下来继续使用 (cuda-gdb) step，看看 __hmul(...) 的实现：

这里的 NV_IF_ELSE_TARGET(cond, , ) 表示可能存在两种可能的实现方式，根绝第一个表达式的真假来选择后边的第二个或者第三个表达式。因为我们使用了 sm_120, 不等于 sm_53，可以初步猜测是调用了后边的第三个表达式的内容来实现乘法。接下来通过 cuda-gdb 来单步调试验证一下。

我们已经猜测会执行后边三行代码：2653，2654，2655等，但是为了验证，这里做了个新函数 hhhaddd()，插入到第三个表达式的中间 float xfa = hhhaddd(fa); ：

这会导致计算结果必然是错误的，但是可以给这个 hhhaddd 打断点，然后直接 continue，果然停在了这个函数上。

说明执行了这三行代码，即，half 的乘法，是使用 float32 的乘法指令来实现的:

    const float fa = __half2float(a);const float fb = __half2float(b);return __float2half(fa * fb);

2. 写个 cpu gemm_fp16

矩阵小一点，方便验证，其中的输出格式，是为了能够简单地放进matlab 做对比验证：


#include <stdio.h>
#include <stdlib.h>
#include "cuda_fp16.h"void init_matrix(half *A, int lda, int m, int n, bool colMajor)
{if(colMajor){for(int j=0; j<n; j++){for(int i=0; i<m; i++){half x = half(rand()*1.0f/RAND_MAX);A[i + j*lda] = x;printf(" %f",  float(x));}}printf("\n\n");}else{for(int i=0; i<m; i++){for(int j=0; j<n; j++){half x = half(rand()*1.0f/RAND_MAX);A[i*lda + j] = x;printf(" %f",  float(x));}}}
}void print_matrix(half *A, int lda, int m, int n, bool colMajor)
{printf("[ ...\n");for(int i=0; i<m; i++){for(int j=0; j<n; j++){if(colMajor)printf(" %5.4f, ", float(A[i + j*lda]));elseprintf(" %5.4f, ", float(A[i*lda + j]));}printf(" ; ...\n");}printf("]\n");
}void gemm_fp16_cpu(int M, int N, int K,half* A, int lda,half* B, int ldb,half* C, int ldc,half alpha, half beta)
{for(int i=0; i<M; i++){for(int j=0; j<N; j++){half sigma = half(0.0);for(int k=0; k<K; k++){sigma += A[i + k*lda]*B[k + j*lda];}C[i + j*ldc] = alpha*sigma + beta*C[i + j*ldc];}}
}int main()
{int m = 4;int n = 4;int k = 4;int lda = m;int ldb = k;int ldc = m;half *A_h;half *B_h;half *C_h;half alpha = half(1.0);half beta  = half(1.0);A_h = (half*)malloc(lda * k * sizeof(half));B_h = (half*)malloc(ldb * n * sizeof(half));C_h = (half*)malloc(ldc * n * sizeof(half));init_matrix(A_h, lda, m, k, true);init_matrix(B_h, ldb, k, n, true);init_matrix(C_h, ldc, m, n, true);printf("A_h =");print_matrix(A_h, lda, m, k, true);printf("B_h =");print_matrix(B_h, ldb, k, n, true);printf("C_h =");print_matrix(C_h, ldc, m, n, true);gemm_fp16_cpu(m, n, k, A_h, lda, B_h, ldb, C_h, ldc, alpha, beta);printf("C_h =");print_matrix(C_h, ldc, m, n, true);return 0;
}

Makefile：

EXE := hello_gemm.fp16all: $(EXE)%: %.cunvcc --gpu-architecture=sm_120 -g $< -o $@ -I /usr/local/cuda/include.PHONY: clean
clean:-rm -rf $(EXE)

编译运行：

$ make

octave 验证：

误差范围内，结果是相等的。

3. GPU 的最简单版本 gemm_v01

简单主要是指没有任何优化考虑。单个warp 工作，也不考虑数据复用、异步加载，不考虑 tensor core 加速，流水线等都不考虑。

我们可以先稍微看看 RTX 5080 的硬件信息：

10752 个cuda core，每个warp 占 32 个 cuda core【注，从 Ampere 开始，每个warp 同时占用 32 个 cuda core；之前架构是 16 个 cuda core 迭代两次完成 32 个 thread 的任务；】，
总共含 84 个sm，
所以，每个sm 存在 128个 cuda core，也就是 128/32 = 4 个同时运行的 warp，也即 4 个 tensor core/sm；也就是每个 block 最多可以同时占用 4 个tensor core。

这个 v01 版本不考虑使用 tensor core，仅启动单个warp 工作。

ex/hello_gemm.fp16.cu


#include <stdio.h>
#include <stdlib.h>
#include "cuda_fp16.h"void init_matrix(half *A, int lda, int m, int n, bool colMajor)
{if(colMajor){for(int j=0; j<n; j++){for(int i=0; i<m; i++){half x = half(rand()*1.0f/RAND_MAX);A[i + j*lda] = x;printf(" %f",  float(x));}}printf("\n\n");}else{for(int i=0; i<m; i++){for(int j=0; j<n; j++){half x = half(rand()*1.0f/RAND_MAX);A[i*lda + j] = x;printf(" %f",  float(x));}}}
}void print_matrix(half *A, int lda, int m, int n, bool colMajor)
{printf("[ ...\n");for(int i=0; i<m; i++){for(int j=0; j<n; j++){if(colMajor)printf(" %5.4f, ", float(A[i + j*lda]));elseprintf(" %5.4f, ", float(A[i*lda + j]));}printf(" ; ...\n");}printf("]\n");
}void gemm_fp16_cpu(int M, int N, int K,half* A, int lda,half* B, int ldb,half* C, int ldc,half alpha, half beta)
{for(int i=0; i<M; i++){for(int j=0; j<N; j++){half sigma = half(0.0);for(int k=0; k<K; k++){sigma += A[i + k*lda]*B[k + j*lda];}C[i + j*ldc] = alpha*sigma + beta*C[i + j*ldc];}}
}__global__ void gemm_v01_fp16_all(int M, int N, int K,half* A, int lda,half* B, int ldb,half* C, int ldc,half alpha, half beta)
{unsigned int i = threadIdx.x;unsigned int j = threadIdx.y;if(i==16) printf("%d ", j);
printf("threadIdx.x=%d  ", threadIdx.x);half sigma = half(0.0);for(unsigned int k = 0; k<K; k++){sigma += A[i + k*lda]*B[k + j*ldb];}C[i + j*ldc] = alpha*sigma + beta*C[i + j*ldc];
}void gemm_v01_test(int m, int n, int k,half* Ah, int lda,half* Bh, int ldb,half* Ch, int ldc,half alpha, half beta,half* Cd2h)
{//1. alloc ABC_dhalf * Ad = nullptr;half * Bd = nullptr;half * Cd = nullptr;cudaMalloc((void**)Ad, lda*k*sizeof(half));cudaMalloc((void**)Bd, ldb*n*sizeof(half));cudaMalloc((void**)Cd, ldc*n*sizeof(half));//2. cpy H2DcudaMemcpy(Ad, Ah, lda*k*sizeof(half), cudaMemcpyHostToDevice);cudaMemcpy(Bd, Bh, ldb*n*sizeof(half), cudaMemcpyHostToDevice);cudaMemcpy(Cd, Ch, ldc*n*sizeof(half), cudaMemcpyHostToDevice);//3. Gemm_v01, simple cuda core gemmdim3 block_;dim3 grid_;block_.x = 32;block_.y = 32;grid_.x = 1;grid_.y = 1;printf("__00________\n");gemm_v01_fp16_all<<<grid_,block_>>>(m, n, k, Ad, lda, Bd, ldb, Cd, ldc, alpha, beta);
printf("##11########\n");//4. cpy D2HcudaMemcpy(Cd2h, Cd, ldc*n*sizeof(half), cudaMemcpyDeviceToHost);//5. free ABC_dcudaFree(Ad);cudaFree(Bd);cudaFree(Cd);
}
int main()
{int m = 32;int n = 32;int k = 32;int lda = m;int ldb = k;int ldc = m;half *A_h;half *B_h;half *C_h;half *C_d2h;half alpha = half(1.0);half beta  = half(1.0);A_h = (half*)malloc(lda * k * sizeof(half));B_h = (half*)malloc(ldb * n * sizeof(half));C_h = (half*)malloc(ldc * n * sizeof(half));C_d2h = (half*)malloc(ldc * n * sizeof(half));init_matrix(A_h, lda, m, k, true);init_matrix(B_h, ldb, k, n, true);init_matrix(C_h, ldc, m, n, true);printf("A_h =");print_matrix(A_h, lda, m, k, true);printf("B_h =");print_matrix(B_h, ldb, k, n, true);printf("C_h =");print_matrix(C_h, ldc, m, n, true);gemm_fp16_cpu(m, n, k, A_h, lda, B_h, ldb, C_h, ldc, alpha, beta);printf("C_h =");print_matrix(C_h, ldc, m, n, true);gemm_v01_test(m, n, k, A_h, lda, B_h, ldb, C_h, ldc, alpha, beta, C_d2h);printf("C_d2h =");print_matrix(C_d2h, ldc, m, n, true);return 0;
}

未完待续

。。。。