• 各位老板好，我今儿要给各位演示一下，在 sm75 机器上用 cuda core计算 fp32 矩阵乘法！
  • 同时和cublas比较性能呢！
• 由于sm75机器上没有tf32 的tensor core，因此在sm75 如T4机器上计算sgemm，只能用cuda core计算！
• 下面我将逐步解释代码，要让小白都能听得懂！

• 估计很多用户想立刻把我写的这个demo先跑起来看看是啥效果，那就请你把我下面的代码放到文件里面
• 文件起个名字叫A.cu
• 然后用 编译命令

nvcc A.cu -o a.out -arch sm_75 -lcublas -I /root/cutlass/include/ -std=c++17
/root/cutlass是cutlass仓库的路径，如果cutlass仓库在别的目录下，请小可爱你自己修改！

代码里面比较了cublas和cutlass的性能，请各位看官自己测试下！我测试的性能是下面这样的，你也快来测试下吧！

T4，fp32，性能

• 下面开始看我写的代代码吧！请用户先看 void init函数以及它之前的部分。
• 首先映入眼帘的是一些头文件和宏定义，其中宏定义WARMUP和REPEATE是用来测试性能的！
• init函数用来初始化矩阵，这里用随机数初始化！

• 紧接着是下面的函数，CutlassSgemmNN，这个函数就是将cutlass提供的device级别的Gemm类以及与其相关的一些函数封装了一下。
• 这个函数假设A，B，C都是column major的哦！请你一定要注意哦！
  • 他的参数有M，N，K alpha，beta，A，B，C，都是好理解的参数
  • lda表示啥呢？
    • （1）如果A是row major，他表示的是A矩阵的第[0][0]和[1][0]之间的物理距离，也就是两行之间同一个元素间的距离！
    • （2）如果A是col major，他表示的是A矩阵的第[0][0]和[0][1]之间的物理距离，也就是两列之间同一个元素间的距离！
      • 有人很好奇，卧槽，这个lda不就应该是M吗？其实不一定哦，你想一下如果A矩阵是某个更大矩阵的子矩阵呢！是不是就不能说他是M啦！
    • 这里由于我们假设A是col major，lda的含义就是(2)了！并且由于我这个例子里面$A$是一个完整的矩阵，因此lda其实就是M啦！

• 下面就是main函数啦！
• 我们构建了一个a和一个b和一个c，我们在cpu上计算$c=a\ast b$（也就是代码里面baseline的计算哦！）时候是将abc都看成row major的矩阵的！
• 可以我们将要调用的CutlassSgemmNN默认输入输出都是col major的呢，这个咋办呀？
• 这里面涉及到一个技巧
  • 那就是row major的A矩阵shape为[M,K]
  • A.T的shape显然是[K,M]
  • 如果A.T也是row major，那么A.T和A的在内存中的线性数据肯定不同的！
  • 如果A.T是col major，那么A.T和A在内存中的线性数据就相同了！
    • 上面这两个数据在内存中线性排布是一摸一样的！！！仔细想想啊！
• 关于矩阵运算在数学上有，如果c=a*b，那么c.T=b.T*a.T
  • 也就是说我们为了求得到row major的c矩阵，可以通过获得col major的c.T来求得！因为他俩在内存中线性排布是一摸一样的！
  • 也就是用col major的b.T和col major的a.T来求得即可！
  • 而col major的 b.T 的地址，显然就是b地址啊！
  • col major的 a.T 的地址，显然就是a地址啊！
• 其他的一些注释我都放到代码里啦！请各位看官看下哦！

#include <stdio.h>
#include <chrono>
#include <ctime>
#include <iostream>
#include <ratio>
#include "cublas_v2.h"
#include "cutlass/gemm/device/gemm.h"#define WARMUP 10
#define REPEATE 10using DATATYPE = float;
using ACCU_DATATYPE = float;void init(DATATYPE *a, int size) {for (int i = 0; i < size; i++) {a[i] = (rand() % 9999) / 10000.0;}
}cudaError_t CutlassSgemmNN(int M, int N, int K, float alpha, float const *A, int lda, float const *B, int ldb, float beta,float *C, int ldc) {// 这个 ColumnMajor 这句话太好理解了，就是col major的意思！using ColumnMajor = cutlass::layout::ColumnMajor;// 下面实例化了一个模版类`class LinearCombination`,这个类的一些参数请看这个文件// https://github.com/NVIDIA/cutlass/blob/main/include/cutlass/epilogue/thread/linear_combination.h吧！using EpilogueOutputOp = cutlass::epilogue::thread::LinearCombination<float, 1, float, float>;// 下面这个牛逼了！实例化了一个模版类`class Gemm`,这个模版类的参数很多，请你们看这个文件吧。// https://github.com/NVIDIA/cutlass/blob/main/include/cutlass/gemm/device/gemm.husing CutlassGemm = cutlass::gemm::device::Gemm<float, ColumnMajor,float, ColumnMajor,float, ColumnMajor,float,cutlass::arch::OpClassSimt,cutlass::arch::Sm75,cutlass::gemm::GemmShape<128,128, 8>,cutlass::gemm::GemmShape<32,64, 8>,cutlass::gemm::GemmShape<1,1, 1>, EpilogueOutputOp,cutlass::gemm::threadblock::GemmIdentityThreadblockSwizzle<>,2,1,1,true >;// 上面只是实例化了模版类，下面这行是真的用模版类实例化了对象！CutlassGemm gemm_operator;// 下面这些是参数哦！// 计算的是D = alpha * A * B + beta * C// 这里面D也就是C哦，又因为beta我们在demo里面给了他是0，alpha给了是1，因此就相当于计算C=A*B啦！CutlassGemm::Arguments args({M, N, K},  // Gemm Problem dimensions{A, lda},   // Tensor-ref for source matrix A{B, ldb},   // Tensor-ref for source matrix B{C, ldc},   // Tensor-ref for source matrix C{C, ldc},   // Tensor-ref for destination matrix D{alpha, beta}, // Scalars used in the Epilogue2); // 这个2表示用split-k算法，且k=2哦！size_t bytes = CutlassGemm::get_workspace_size(args);void * workspace;cudaMalloc( (void**)&workspace, bytes);cutlass::Status status = gemm_operator(args, workspace);if (status != cutlass::Status::kSuccess) {return cudaErrorUnknown;}return cudaSuccess;
}int main(void) {int m = 512;int n = 512;int k = 512;DATATYPE *a, *b;a = (DATATYPE *)malloc(sizeof(DATATYPE) * m * k);b = (DATATYPE *)malloc(sizeof(DATATYPE) * k * n);init(a, m * k);init(b, k * n);ACCU_DATATYPE *c;c = (ACCU_DATATYPE *)malloc(sizeof(ACCU_DATATYPE) * m * n);memset(c, 0, sizeof(ACCU_DATATYPE) * m * n);float *c_cpu_fp32 = (float *)malloc(sizeof(float) * m * n);memset(c_cpu_fp32, 0, sizeof(float) * m * n);DATATYPE *dev_a, *dev_b;ACCU_DATATYPE *dev_c;cublasHandle_t handle;cublasCreate(&handle);// allocate the memory on the GPU cudaMalloc((void **)&dev_a, m * k * sizeof(DATATYPE));cudaMalloc((void **)&dev_b, k * n * sizeof(DATATYPE));cudaMalloc((void **)&dev_c, m * n * sizeof(ACCU_DATATYPE));cudaMemcpy(dev_a, a, m * k * sizeof(DATATYPE), cudaMemcpyHostToDevice);cudaMemcpy(dev_b, b, k * n * sizeof(DATATYPE), cudaMemcpyHostToDevice);cudaEvent_t beg, end;for (int i = 0; i < WARMUP + REPEATE; i++) {if (i == WARMUP) {cudaEventCreate(&beg);cudaEventCreate(&end);cudaEventRecord(beg);}const float alpha = 1.0f;const float beta = 0.0f;// 这里为了求得`row major`的`c`矩阵，选择计算`col major`的`c.T`，因为这俩个地址是一模一样的！// `c.T=b.T*a.T`我们为了获得col major的c.T，// `col major`的`b.T` 的shape是[n,k]，地址就是dev_b// `col major`的`a.T` 的shape是[k,m],地址就是dev_a// 因此下面的参数是下面填写的那样哦！cublasSgemm(handle,CUBLAS_OP_N,CUBLAS_OP_N,n,m,k,&alpha,dev_b,n,dev_a,k,&beta,dev_c,n);// 这个是cutlass的gemm！// CutlassSgemmNN(n, m, k, alpha, dev_b, n, dev_a, k, beta, dev_c, n);}cudaEventRecord(end);cudaEventSynchronize(end);float elapsed_time;cudaEventElapsedTime(&elapsed_time, beg, end);printf("gpu gemm compute time: %f ms\n", elapsed_time);// 把gpu运算的结果拷贝到host端c上！cudaMemcpy(c, dev_c, m * n * sizeof(ACCU_DATATYPE), cudaMemcpyDeviceToHost);// 在cpu上计算结果，这个结果作为baseline，用来确保cutlass没有算错哦！for (int i = 0; i < m; i++) {for (int j = 0; j < n; j++) {double sum = 0.f;for (int ii = 0; ii < k; ii++) {sum += a[i * k + ii] * b[ii * n + j];}c_cpu_fp32[i * n + j] = sum;}}// 看看baseline和gpu上的diff，然后输出！double max_diff = -1.;for (int i = 0; i < m; i++) {for (int j = 0; j < n; j++) {double c_gpu_fp32 = c[i * n + j];if (std::abs(c_cpu_fp32[i * n + j] - c_gpu_fp32) > max_diff) {max_diff = std::abs(c_cpu_fp32[i * n + j] - c_gpu_fp32);}}}printf("max_diff: %f\n", max_diff);cudaDeviceReset();free(a);free(b);free(c);free(c_cpu_fp32);return 0;
}