简介
NVIDIA 率先在 GPU 中引入了通用计算能力，使得开发者能利用 CUDA 编程语言来驱动。这时候 GPU 的核心都是 CUDA Core。由于一个 GPU 里面有大量的 CUDA Core，使得并行度高的程序获得了极大的并行加速。但是，CUDA Core 在一个时钟周期只能完成一个操作，矩阵乘法操作依然需要耗费大量的时间。NVIDIA 为了进一步加速“加乘运算”，在 2017 年推出了 Volta 架构的 GPU，从这个架构开始 Tensor Core 被引入。它可以在一个时钟周期完成两个 4×4x4 半精度浮点矩阵的乘法（64 GEMM per clock）。
基本概念
CUDA编程指北：从入门到实践 未读完，可继续。 CUDA 程序一般使用 .cu 后缀，编译 CUDA 程序则使用 nvcc 编译器。一般而言，一个 CUDA 程序的结构如下：
int main() {
主机代码; // 负责 CPU 和 GPU数据传输、GPU内存管理、以及启动 GPU 内核（内核启动参数指定了 GPU 上线程的数量和分布方式）等
核函数调用; // 每个内核函数在 GPU 的众多 CUDA 核心上并行执行，在 GPU 的多个线程上同时执行
主机代码;
核函数调用;
......
return 0;
}
__global__ void 核函数1(parameters) {
// 在设备代码中，计算任务被分解为多个线程，这些线程组成线程块（Block），多个线程块组成一个线程网格（Grid）。CUDA 提供了 threadIdx、blockIdx 等内置变量来获取线程的索引，从而让每个线程在数据中找到属于自己的计算任务。
......
}
__global__ void 核函数2(parameters) {
......
}
前缀__global__用来定义一个核函数，在 CUDA 中，核函数只能返回 void 类型（无返回值），这意味着当我们需要写计算结果时，应该在参数列表中传入一个用来存放计算结果的指针，然后将计算结果写回到这个指针指向的存储空间中。CUDA 核函数传入的参数必须是指向设备内存，因此，我们必须预先在主机代码中分配设备内存并初始化。分配设备内存可以使用 cudaMalloc 函数，初始化设备内存则可以将一段已经初始化好的主机内存拷贝到刚分配好的设备内存中，这可以使用 cudaMemcpy 函数实现，这两个函数的函数原型如下：
cudaError_t cudaMalloc(void** d_ptr, unsigned int size);
cudaError_t cudaMemcpy(void* d_ptr, void* h_ptr, unsigned int size, enum cudaMemcpyKind)
PS：所以在推理框架中，显存管理是推理框架负责。核函数 都是被封装后，注册到算子里，被类似op.forward 触发执行。 核函数是无状态的。 核函数的调用语法（内核启动语法）如下所示：
// CUDA 使用特殊的语法 <<<Grid, Block>>> 启动内核函数。
kernel_function<<<grid_size, block_size>>>(parameters)
// 也可以认为是
kernel<<<numBlocks, threadsPerBlock>>>(parameters)
// numBlocks 表示线程块的数量，threadsPerBlock 表示每个线程块中包含的线程数。
// 通过指定线程块数和线程数，内核启动控制了 GPU 的并行粒度。较大的数据集通常需要更多的线程和线程块来充分利用 GPU 的并行能力。
CUDA 的核函数设计一般遵循如下范式：data1，data2 … 表示需要处理的数据指针，index1 和 index2 … 用来定位需要计算的数据的位置，some_operation 对这些数据进行指定的计算操作，然后写回到参数列表中传入的用于记录结果的 result 指针中。总结下来就是两部曲：确定线程和数据的对应； 对需要处理的数据执行操作。PS：CUDA 最难的是并行思想。并行思想其中难上加难的东西是数据分组。并行计算中，最重要的一点是为数据分组成多个小数据块，每个线程（进程）再去实现SPMD或者SIMD/T。而这个数据分组的方式，存储方法等等直接的影响到你这个并行程序最终的性能。大部分的并行程序，解决了数据分组问题，其本身的问题就解决了，算法本身的优化反倒是不是那么的重要了。
__global__ void kernel_function(data1, data2, ..., result) {
index1, index2, ... = get_index(thread_info)
result = some_operations(data1[index1], data2[index2], ...)
}
内核启动后，GPU 可以异步执行任务，CPU 继续进行其他操作，直至需要等待 GPU 完成。开发者可以利用这种异步特性，使程序在 CPU 和 GPU 间并行执行，达到更高的并行效率。此外，CUDA 提供了同步函数（如 cudaDeviceSynchronize），确保 CPU 在需要时等待 GPU 完成所有操作，避免数据不一致的问题。
编译 CUDA 程序：编译 CUDA 程序需要使用 Nvidia 官方提供的编译器 nvcc。nvcc 会先将所有源代码先分离成主机代码和设备代码，主机代码完整支持 C++ 语法，设备代码只部分支持 C++ 语法。nvcc 先将设备代码编译为 PTX（parallel thread execution）伪汇编代码，再将 PTX 代码编译为二进制的 cubin 目标代码。CUDA 中核函数也因此不能直接作为类的成员函数，如果希望使用面向对象，我们一般通过包装函数调用核函数，然后将这个包装函数作为成员函数。
cuda的寄存器存储和共享内存
// 静态声明shared memory
__global__ void my_kernel() {
__shared__ int i;
}
// 动态声明shared memory
my_kernel(grid_dim, block_dim, 8) // 调用的时候指定动态内存的大小。这里8就是动态分配的。
__global__ void my_kernel() {
extern __shared__ int arr[];
}
// kernel函数中寄存器分配的场景一般有如下几种：临时变量，循环中分配的变量，函数调用和函数返回值。
__device__ int do_something(int x) {
return x*100; //返回值会被存放在register上。
}
__global__ void my_kernel() {
int a; //临时变量.a会被存放在register中
// for中的变量i会被存放在register中
for(int i=0;i<100;++i){
...
}
// do_something()的返回值会被存放在register中
int b = do_something(a);
// 大数组不会被存放在register上
int arr[1000];
}
大数组、结构体、动态分配的数组不会被存放在register上，而是存放在local memory中。local memory是global memory中的一块区域，由cuda自动分配，专门用来存放线程私有的数据。很显然，它的访问速度会比register 和shared memory慢很多。
cuda的同步有两种：
- 系统级同步：同步host和device的工作，用cudaDeviceSynchronize()接口。这个接口会阻塞所有host的工作，直到cuda端的工作完成。
- block线程同步。同步同一个block内的线程，用__synthreads()接口。同一个block内的线程用register和shared memory 进行通信。 cuda不同block之间的线程无法同步。如果需要，只能使用系统级同步方式，，使用cudaDeviceSynchronize()进行等待，在不同block的线程达到checkpoint后结束当前的kernel，开启新的kernel。
执行过程
CUDA C++ 编程指北-第一章：入门以及编程模型 未读。
device 函数和global函数因为需要在GPU上运行，因此不能调用常见的一些 C/C++ 函数（因为这些函数没有对应的 GPU 实现）。
| 限定符 | 执行 | 调用 | 备注 |
|---|---|---|---|
| global | 设备端执行 | 可以从主机调用也可以从某些特定设备调用 | 异步操作，host 将并行计算任务发射到GPU的任务调用单之后，不会等待kernel执行完就执行下一步 |
| device | 设备端执行 | 设备端调用 | |
| host | 主机端执行 | 主机调用 |
典型的CUDA程序的执行流程如下：
- 分配host内存，并进行数据初始化；
- 分配device内存，并从host将数据拷贝到device上；
- 调用CUDA的核函数在device上完成指定的运算；
- 将device上的运算结果拷贝到host上；
- 释放device和host上分配的内存。
矩阵加法示例
// __global__ 表示在device上执行，从host中调用
// 两个向量加法kernel，grid和block均为一维
// 为每个线程分配了一个其所属的向量元素，然后驱动线程分别完成计算
__global__ void add(float* x, float * y, float* z, int n){
// 计算线程的编号，根据自己的线程编号，读取数据源不同位置的元素，执行计算，并将结果写入结果的不同位置。
// 这样，我们就为不同线程安排了独立的工作，让他们并发地完成工作。
int index = threadIdx.x + blockIdx.x * blockDim.x;
// 步长
int stride = blockDim.x * gridDim.x;
for (int i = index; i < n; i += stride){ // cuda编程经常会跨步长，每个thread 只计算数据的局部和/积/xx。
z[i] = x[i] + y[i]; // 每个线程都执行这条指令，每个线程读取不同元素执行相同计算
}
}
int main(){
int N = 1 << 20;
int nBytes = N * sizeof(float);
// 申请host内存
float *x, *y, *z;
x = (float*)malloc(nBytes);
y = (float*)malloc(nBytes);
z = (float*)malloc(nBytes);
// 初始化数据
for (int i = 0; i < N; ++i){
x[i] = 10.0;
y[i] = 20.0;
}
// 申请device内存
float *d_x, *d_y, *d_z;
cudaMalloc((void**)&d_x, nBytes);
cudaMalloc((void**)&d_y, nBytes);
cudaMalloc((void**)&d_z, nBytes);
// 将host数据拷贝到device
cudaMemcpy((void*)d_x, (void*)x, nBytes, cudaMemcpyHostToDevice);
cudaMemcpy((void*)d_y, (void*)y, nBytes, cudaMemcpyHostToDevice);
// 定义kernel的执行配置
dim3 blockSize(256);
dim3 gridSize((N + blockSize.x - 1) / blockSize.x);
// 执行kernel
add << < gridSize, blockSize >> >(d_x, d_y, d_z, N); # 第一个数字指明改程序分配多少个block，第二个数字程序指明每个block中的thread个数
// 将device得到的结果拷贝到host
cudaMemcpy((void*)z, (void*)d_z, nBytes, cudaMemcpyDeviceToHost);
// 检查执行结果
float maxError = 0.0;
for (int i = 0; i < N; i++)
maxError = fmax(maxError, fabs(z[i] - 30.0));
std::cout << "最大误差: " << maxError << std::endl;
// 释放device内存
cudaFree(d_x);
cudaFree(d_y);
cudaFree(d_z);
// 释放host内存
free(x);
free(y);
free(z);
return 0;
}
如何在 CPU 之上调用 GPU 操作？可以通过调用 __global__ 方法来在GPU之上执行并行操作。我的第一份CUDA代码 - xcyuyuyu的文章 - 知乎 kernel在调用时也必须通过执行配置<<<grid, block>>>来指定kernel所使用的线程数及结构。一个thread需要两个内置的坐标变量（blockIdx，threadIdx）来唯一标识，其中blockIdx指明block在grid中的位置，而threaIdx指明线程所在block中的位置。
main函数中的cudaMalloc、cudaMemcpy，是CPU操作GPU内存的操作，在分离式GPU架构（也就是独显）中，CPU分配内存用于GPU计算，再将数据传输到分配的内存空间，然后在GPU上启动内核函数。GPU执行的内核函数只能从分配的GPU内存空间读取数据。代码中的host向量对应CPU内存的数据，而device向量则代表GPU内存的数据。完成内存分配和数据拷贝后，CPU触发GPU执行add内核函数。触发时同时指定了执行内核函数的线程的组织形式。在CUDA编程中，线程以thread，thread block，grid的层级结构进行组织。一个线程块由多少个线程组成可以指定，线程块本身的数量则是由计算规模决定的，比如根据向量的长度计算了线程块的数量：int blocks_num = (n + t - 1) / t; // thread block数量，这样计算的目的是保证线程数量足够，即每一个计算单元都有一个线程负责计算。add内核函数就是以SIMT模型进行编程的，安排所有线程执行相同的指令，但每个线程执行指令时的指令操作数均不同，这便是SIMT。
指令集与编译
上面CUDA C语言编写的add 只是到了高级语言层面，众所周知，高级语言需要转换为机器码才能被机器执行，本节将简单介绍CUDA C/C++的程序的编译流程，以及CUDA的PTX、SASS指令集。
- SASS（Streaming Assembly）是GPU的机器指令集，是实际在GPU上执行的指令。SASS指令集直接对应GPU架构（Maxwell、Pascal等），虽然不是严格的一一对应，但通常每个GPU架构有专属的SASS指令集，因此需要针对特定架构进行编译。
- PTX（Parallel Thread Execution）是一种中间表示形式，位于高级GPU编程语言（如CUDA C/C++）和低级机器指令集（SASS）之间。PTX与GPU架构基本无耦合关系，它本质上是从SASS上抽象出来的一种更上层的软件编程模型，PTX的存在保证了代码的可移植性（同一份PTX分发到不同架构上转为对应SASS）与向后兼容性（可将PTX代码转为最新GPU架构对应的SASS）。PTX是开发者可编程的最底层级，而SASS层则是完全闭源的，这也是NVIDIA的“护城河”之一。
CUDA程序的编译由NVCC（NVIDIA CUDA Compiler）完成。
首先，NVCC完成预处理；随后分类代码为设备代码和主机代码，NVCC驱动传统的C/C++编译器主机代码的编译和汇编；对于设备代码，NVCC将其编译针对某架构的SASS，编译过程中涉及C –> PTX –> SASS的转化，但通常不显式表现出来，生成的PTX/SASS码也会被直接嵌入最终的可执行文件。
运行期，GPU会优先查找可执行文件中是否有适合当前架构的SASS，如有则直接执行。若无，则GPU驱动（driver）会使用JIT（Just-In-Time）编译手段，将PTX码编译为当前架构对应的SASS再执行（前提是可执行文件必须包含PTX）。
矩阵计算优化（没太懂）
GPU存储可分为物理内存（硬件真实存在的）和逻辑内存（由cuda做抽象的）。 为什么要这么分呢？因为各个GPU的物理内存架构是不一样的，如果你写代码时还要考虑每个GPU的独特性，那可太痛苦了。所以cuda在这里帮了大忙：它对内存架构做了一层抽象，你只要按照它抽象后的框架写代码就可以。
每个thread占用一个SP（cuda core），即1个warp会占用1个SM上的32个SP。 有了这些前置知识，现在我们可以来看cuda矩阵优化的过程了。 假设矩阵
- A = (M,K) = (512,512)
- B = (K,N) = (512,512)
- C = AB = (M,K) * (K,N) = (512,512)
每个thread负责读取A矩阵的一行和B矩阵的一列，去计算C矩阵的一个元素。则一共需要M*N个thread。 矩阵A和矩阵B都存储在global memory，每个thread直接从global memory上进行读数，完成计算：
- 为了计算出C中的某个元素，每个thread每次都需要从global memory上读取A矩阵的一行（K个元素），B矩阵的一列（K个元素），则每个thread从global memory上的读取次数为2K
- C中共有M*N个thread，则为了计算出C，对global memory的总读取次数为： 2MNK
Naive GEMM的代码见下（完整代码见 sgemm_naive.cu ）：
// 将二维数组的行列索引转成一维数组的行列索引，这样可以更高效访问数据
// row, col：二维数组实际的行列索引，ld表示该数组实际的列数
// 例：二维数组实际的行列索引为(1, 3)，即第二行第四个元素，二维数据的总列数 = 5
// 返回的一位数组形式的索引为: 1*5 + 3 = 8
#define OFFSET(row, col, ld) ((row) * (ld) + (col))
// 定义naive gemm的kernel函数
__global__ void naiveSgemm(float * __restrict__ a, float * __restrict__ b, float * __restrict__ c, const int M, const int N, const int K) {
// 当前thread在C矩阵中的row
int m = blockIdx.y * blockDim.y + threadIdx.y;
// 当前thread在C矩阵中的col
int n = blockIdx.x * blockDim.x + threadIdx.x;
if (m < M && n < N) {
float psum = 0.0;
// 告知编译器自动展开循环体，这样可以减少循环控制的开销（循环次数小的时候可以这么做）
#pragma unroll
// 取出A[row]和B[col]，然后逐个元素相乘累加，得到最终结果
for (int k = 0; k < K; k++) {
// a[OFFSET(m, k, K)]: 获取A[m][k]
// b[OFFSET(k, n, N)]: 获取B[k][n]
psum += a[OFFSET(m, k, K)] * b[OFFSET(k, n, N)];
}
c[OFFSET(m, n, N)] = psum;
}
}
const int BM = 32, BN = 32;
const int M = 512, N = 512, K = 512;
dim3 blockDim(BN, BM);
dim3 gridDim((N + BN - 1) / BN,(M + BM - 1) / BM);
GPU/CUDA/驱动和机器学习训练框架的关系
显卡是硬件，硬件需要驱动，否则不能调用其计算资源。CUDA又是什么？
- 在2007年之前，GPU由CPU操作，CPU把一些图形图像的计算任务交给GPU执行。程序员不需要与GPU打交道。随着GPU计算能力的发展，越来越多的计算场景由GPU完成效果会更好。但现有的程序无法直接自由控制GPU的处理器。当然程序员也可以直接写代码与显卡驱动对接，从而直接控制GPU的处理器，但这样代码恐怕写起来要让人疯掉。nvidia当然会有动力提供一套软件接口来简化操作GPU的处理器。nvidia把这一套软件定义为CUDA。
- 多核 CPU 和众核 GPU 的出现意味着主流处理器芯片现在是并行系统。挑战在于开发能够透明地扩展可并行的应用软件，来利用不断增加的处理器内核数量。CUDA 并行编程模型旨在克服这一挑战，同时为熟悉 C 等标准编程语言的程序员保持较低的学习曲线。CUDA 编程手册系列第一章：CUDA 简介
gpu 和 cuda 和 gpu driver 之间的关系：比如 TX3090需要Compute Capability在8.6以上的cuda，而满足这个要求的cuda又只有11.0以上的版本。而cuda11版本又需要版本号>450的显卡驱动。
显卡
==> Compute Capability 查看显卡支持的Compute Capability, https://developer.nvidia.com/cuda-gpus
==> cuda
==> GPU driver 查看cuda对驱动的要求 (https://docs.nvidia.com/cuda/cuda-toolkit-release-notes/index.html
==> tf/pytorch version
Compute Capability的数值和GPU的计算速度无关，但是和GPU可执行的任务种类有关。The compute capability of a device is represented by a version number, also sometimes called its “SM version”. This version number identifies the features supported by the GPU hardware and is used by applications at runtime to determine which hardware features and/or instructions are available on the present GPU.
从模型到算子，以卷积计算为例
- 从卷积到矩阵乘 将输入数据（Feature Map）和卷积核数据进行重排，卷积操作本质上可以等效理解为矩阵乘操作。卷积操作的过程大概可以描述为按照约定的窗口大小和步长，在 Feature Map 上进行不断地滑动取数，窗口内的 Feature Map 和卷积核进行逐元素相乘，再把相乘的结果累加求和得到输出 Feature Map 的每个元素结果。
- 矩阵乘分块 Tilling。卷积转换后的矩阵乘的维度非常大，而芯片里的内存空间往往是有限的（成本高），表现为越靠近计算单元，带宽越快，内存越小。为了平衡计算和内存加载的时间，让算力利用率最大化，AI 芯片往往会进行由远到近，多级内存层级的设计方式，达到数据复用和空间换时间的效果。根据这样的设计，矩阵乘实际的数据加载和计算过程将进行分块 Tilling 处理。Tiling（平铺）是一种优化技术，它涉及将大的矩阵分解成更小的块或“瓦片”（tiles），这些小块的大小通常与CPU或GPU的缓存大小相匹配，以便可以完全加载到缓存中。
- 矩阵乘的库。矩阵乘作为 AI 模型中的重要性能算子，CPU 和 GPU 的平台上都有专门对其进行优化实现的库函数。比如 CPU 的 OpenBLAS, Intel MKL 等，GPU 的 cuBLAS, cuDNN 等。实现的方法主要有 Loop 循环优化 (Loop Tiling)和多级缓存 (Memory Hierarchy)。
- 矩阵乘的优化。在具体的 AI 芯片或其它专用芯片里面，对矩阵乘的优化实现主要就是减少指令开销，可以表现为两个方面：
- 让每个指令执行更多的 MACs 计算。比如 CPU 上的 SIMD/Vector 指令，GPU 上的 SIMT/Tensor 指令，NPU 上 SIMD/Tensor,Vector 指令的设计。
- 在不增加内存带宽的前提下，单时钟周期内执行更多的 MACs。比如英伟达的 Tensor Core 中支持低比特计算的设计，对每个 cycle 执行 512bit 数据的带宽前提下，可以执行 64 个 8bit 的 MACs，大于执行 16 个 32bit 的 MACs。