GPU Computing

background

GPU Computing：the principal motivation for using the GPU is the prospect of high performance at a relatively low cost
CUDA：CUDA provides a flexible programming model and C-like language for implementing data-parallel algorithms on the GPU
FLOPs：每秒所执行的浮点运算次数（floating-point operations per second）
- https://www.zhihu.com/zvideo/1421421171497304064?playTime=463.3
- compute intensity
- memory latency
  - 我们的GNPU，0.75V/900M，1个时钟cycle近似于1ns，所以一个10w cycle的softmax的计算latency是1ms
  - chip size：做小是为了降低电流传输距离，实现更短的latency
  - 带宽memory bandwidth：例如传输带宽是131G/sec，clock latency是89ns，可以算出每个cycle最多可以传输11659bytes的数据，如果我们的程序只要求了f64的x和f64的y也就是8 bytes的数据量，那么内存利用率就是0.14%
- threads
  - 把任务拆解成高并发，比如可以将矩阵乘法拆解成
  - GPU的设计核心就是支持超多线程，而且线程切换很快
- reg & cache
  - GPU通过设计很多寄存器的方式，维持高效的数据传输效率
  - 不同层次的数据存储，reg/cache/shared mem/global mem
- Cuda Core & Tensor Core
  - 都是基本计算单元
  - Cuda Core算浮点，Tensor Core算int
CUDA kernel：https://www.youtube.com/watch?v=jYCxVirq4d0
- cuda kernel的code写的是1个thread要做的事
- 但是它manipulate两个层次，1是thread block，thread block内部的threads has shared memory，2是thread block构成的launching grid，对应地，GPU的memory hierarchy也分为两级，每个thread的私有mem叫local registers，每个block的mem叫shared memory，除此以外tensor.cuda()放的地方叫global memory
- a kernel executes a grid
- threads使用shared memory的时候要sync：所有thread finish以后整个block才算finish
- parallel的程度取决于具体算子：
  - 如果c[i]=f(c[i-1])，比如cumadd，那这种就不好并行化了
  - 如果是ele-wise的，那就可以高度并行化
1
2
3
4
5
6
// implement cuda kernel of C=A+B
__global__ void VecAddKernel(float* A, float* B, float* C, int n) {
int i = blockDim.x * blockIdx.x + threadIdx.x;
if (i<n)
C[i] = A[i] + B[i]
}

CPU host

在cpu上要运行code，to launch the CUDA kernel
指定/计算threads和blocks数量

void VecAddCPU(float* A, float* B, float* C, int n) {
    int threads_per_block = 512;
    int nblocks = ceil(n/threads_per_block);
    VecAddKernel<<<nblocks, threads_per_block>>>(A, B, C, n);
}

CPU & GPU

#include <stdio.h>

// 在cpu上执行
int main(void)
{
    printf("Hello World from CPU\n");
    return 0;
}


// 在gpu上执行
__global__ void HelloFromGPU(void)
{
    printf("Hello from GPU\n");
}
// 在cpu上发起调用
int main(void)
{
    printf("Hello from CPU\n");
    HelloFromGPU<<<1, 5>>>();
    cudaDeviceReset();
    return 0;
}

__global__：用它去修饰函数HelloFromGPU，成为内核函数，该函数在host端异步调用，在device端执行，无返回值
__device__：限定函数在device上调用，在device上执行，用于gpu内部调用
__host__：限定函数在host上调用，在host上执行，即默认cpu执行，所以这个修饰可以省略
<<>>：也指明了是从host端到device端的内核函数调用，里面的参数是执行配置，例子中kernel函数将被执行5个线程执行

流式多处理器(Streaming Multiprocessors，SM)
- GPU架构是通过复制这种架构的构建来实现硬件并行
  - 1个GPU包含多个SM，每个SM支持许多线程并发执行
  - CUDA采用单指令多线程(Single-Instruction Multiple-Thread，SIMT)来管理和执行GPU上的众多线程
    - 两级线程结构
    - 1个内核的所有线程为一个线程grid，所有线程共享全局内存空间
    - 1个grid由多个线程块（block）组成，一个线程块包含一组线程，同一线程块内的线程（thread）通过同步和共享内存的方式实现协作
    - 不同块内的线程不能协作
  - 当host发布一个kernel时，这个内核网格的线程块就被分配到可用的SM上来执行，多个线程在SM上并发执行，多个线程块可以并发地在一个SM上执行，当线程块终止时，新的线程块又可以在腾出的SM上启动执行
- SM核心组件
  - 核心
  - 共享内存/一级缓存
  - 寄存器文件
  - 加载/存储单元
  - 特殊功能单元
  - 线程束调度器warp
- 并行
  - 任务并行：将任务切分成子任务，多线执行，网格并发是任务并行
  - 数据并行：每个线程处理一份数据，线程执行的计算任务是一样的，线程并行是数据并行
- 线程束
  - CUDA采用SIMT架构来管理和执行线程，将线程块中的线程每32个为一组进行划分，每一组被称为一个线程束(warp)
  - 线程束的大小warpSize是CUDA中的一个内部属性
  - 线程束是GPU的基本执行单元
  - 一旦线程块被调度到一个SM上，线程块中的线程就会被进一步划分为线程束，每个线程束中的所有线程执行相同的命令，每个线程拥有自己的指令地址计数器和寄存器状态
CUDA中的可编程内存的类型
- overview
  - 寄存器(Registers)
  - 本地内存(Local Memory)
  - 共享内存(Shared Memory)
  - 常量内存(Constant Memory)
  - 纹理内存(Texture Memory)
  - 全局内存(Global Memory)
- thread独享：寄存器(Registers) & 本地内存(Local Memory)
- 一个block中的thread共享：共享内存(Shared Memory)
- 其余全局共享
- 寄存器
  - 内核函数中声明且没有其他修饰符修饰的变量，如VectorAddGPU中的线程索引变量i
  - 用于存放内核函数中需要频繁访问的线程私有变量
  - 变量与内核函数的生命周期相同
  - 寄存器是GPU中访问速度最快的内存空间
  - 寄存器溢出：一旦内核函数使用了超过硬件限制的寄存器数量，则会使用本地内存来代替多占用的寄存器
- 共享内存
  - 被__shared__修饰符修饰的变量被存储到共享内存中
  - 生命周期伴随整个线程块
  - 线程通过使用共享内存中的数据可以实现互相之间的协作
  - 使用共享内存必须调用如下函数进行同步：void __sybcthreads()
- 常量内存
  - 常量变量用__constant__修饰符进行修饰
  - 必须在全局空间内和所有内核函数之外进行声明
  - 用于线程束中的所有线程都需要从相同的内存地址中读取数据的情况
- 全局内存
  - GPU中容量最大、延迟最高的内存空间，其作用域和生命空间都是全局的
  - 可以在任何SM设备中被访问到
GPU缓存
- 缓存：
  - 缓存是cpu和内存之间，内存中访问频率高的数据的复制品（映射）
  - CPU找数据或指令的顺序是：先到一级缓存中找，找不到再到二级缓存中找，还找不到就到内存中找
- overview
  - 一级缓存
  - 二级缓存
  - 只读常量缓存
  - 只读纹理缓存
- 每个SM都有一个一级缓存，所有SM共享一个二级缓存
- 一级和二级缓存用来存储本地内存和全局内存中的数据，包括寄存器溢出的部分
Host内存：https://mp.weixin.qq.com/s?__biz=MzkwNjE2ODMyMQ==&mid=2247484370&idx=1&sn=a23ddab4c1f7bff4f221f7ae78048d0f&chksm=c0edd2bcf79a5baa9889b06e721f2c70d62bcc489490e8bb3fdc3a74e4ad77803c447351645f&scene=178&cur_album_id=2132879019699240961#rd
- 可分页内存pageable：指令/数据被离散存储至内存中，因为cpu和gpu是异步的，可能会发生cpu在gpu执行过程中修改内存发生不安全的情况，所以cuda会将临时分配页面锁定
- 固定内存pinned：cuda可以分配固定的主机内存（cudaMallocHost），不需要中间freeze的过程，直接被device端访问，因此可以用很高的带宽进行读写操作
CUDA Stream
- CUDA流是一系列异步操作的集合，同一个CUDA流中的操作严格按照顺序在GPU上运行
- 多个流同时启动多个内核任务就可以实现网格级并发，任务并行
- 显式地创建多个流，去执行数据拷贝、kernel计算等，不同的CUDA操作是可以重叠进行的
- 创建显式流：
  1
  2
  cudaStream_t stream;
  cudaStreamCreate(&stream);
- 销毁流：
  1
  cudaError_t cudaStreamDestroy(cudaStream_t stream);
- 异步拷贝数据：
  1
  cudaError_t cudaMemcpyAsync(void* dst, const void* src, size_t count, cudaMemcpyKind kind, cudaStream_t stream = 0);
- 执行指定流的kernel
  1
  kernel_name<<<grid, block, sharedMemSize, stream>>>(...);
- 考虑矩阵乘法
  - 假如任务是两个矩阵乘法，可以将每个乘法任务绑定一个stream
  - 每个stream异步计算，然后结果同步
tensorRT
- NVIDIA TensorRT是一种高性能神经网络推理(Inference)引擎，用于在生产环境中部署深度学习应用，可提供最大的推理吞吐量和效率
- TensorRT有如下功能
  - Quantization
    - QAT：Quantization Aware Training，训练中量化，通过在浮点训练过程中插入虚拟的int8量化节点
    - PTQ：Post-Training Quantization，训练后量化，那一些样本图片分析不同激活层的结果分布，在以有模型上添加scale实现量化，tensorRT的PTQ以速度优先，无法控制某一层的精度
  - Kernel Auto Tuning：比如一个卷积有多种实现方案，tensorRT会根据实际的部署设备、tensor尺寸、batchsize等参数选择最efficient的方案
  - Elimination of Redundant Layers and Operations：
    - 读写tensor有个memory bandwidth bottleneck，所以layer fusion可以加速
    - 像concat的这种可以预留一个target output buffers直接把相关结果写进去，这是通过layer elimination来实现加速
    - 还有一种情况是shared structure but different weights，tensorRT会把它合并成一个宽的层，类似group convolution的原理，一次层推理得到3个path的结果，这也是layer fusion来使用GPU的计算能力
  - Dynamic Tensor Memory
- 使用流程
  - given a trained model
  - use TensorRT to parse a trained model，做如上优化得到一个plan file，a plan file不仅包含序列化的图，还包含了inference schedule，用什么kernel、什么执行顺序之类的，可以理解为tensorRT Graph
  - load and deserialize a saved plan file to create a TensorRT engine object and run inferenc
- TensorRT配套周边：https://github.com/NVIDIA/TensorRT/tree/master/tools
  - ONNX GraphSurgeon：修改ONNX模型，增加或者剪掉某些节点，修改名字或者维度等等
  - PyTorch-Quantization：在Pytorch训练或者推理的时候加入模拟量化操作，支持量化训练后的模型导出ONNX和TRT
- cookbook：https://github.com/NVIDIA/trt-samples-for-hackathon-cn/tree/master/cookbook
- workflow
  - 使用深度学习框架自带的tensorRT接口，遇到不支持的算子：返回原框架计算
  - 【推荐】导出onnx，然后导入tensorRT，遇到不支持的算子：custom node
  - 使用tensorRT搭建网络，遇到不支持的算子：写plugin(cudac++)
- 命令行工具
  - trtexec：由onnx构建模型序列plan，推理plan及查看相关信息
  - polygraphy：
Jetson
- NVIDIA将TensorRT、cuDNN、CUDA等功能整合到Jetson平台，通过一体化软件包JetPack按需提供给开发者们

pytorch custom cuda op

custom add：https://mp.weixin.qq.com/s?__biz=MzkwNjE2ODMyMQ==&mid=2247484261&idx=1&sn=d7b500fb8e81da96e4828b36d64a749f&chksm=c0edd20bf79a5b1d54e6557be02552c8ab16d7225202bb87e3a89bbd203d0b3fcd643f4be6a6&scene=178&cur_album_id=2132879019699240961#rd

cpu code

#include <iostream>

void VectorAddCPU(const float *const a, const float *const b, float *const c,
                  const int n) {
  for (int i = 0; i < n; ++i) {
    c[i] = a[i] + b[i];
  }
}

int main(void) {
  // alloc memory for host
  const size_t size = 1024;
  float *ha = new float[size]();
  float *hb = new float[size]();
  float *hc = new float[size]();
  for (int i = 0; i < size; ++i) {
    ha[i] = i;
    hb[i] = size - i;
  }

  VectorAddCPU(ha, hb, hc, size);

  delete[] ha;
  delete[] hb;
  delete[] hc;

  return 0;
}

gpu code

#include <cuda_runtime.h>
#include <iostream>

__global__ void VectorAddGPU(const float *const a, const float *const b,
                             float *const c, const int n) {
  int i = blockDim.x * blockIdx.x + threadIdx.x;
  if (i < n) {
    c[i] = a[i] + b[i];
  }
}

int main(void) {
  // 分配CPU内存
  const size_t size = 1024;
  float *ha = new float[size]();
  float *hb = new float[size]();
  float *hc = new float[size]();
  for (int i = 0; i < size; ++i) {
    ha[i] = i;
    hb[i] = size - i;
  }

  // 分配GPU内存
  float *da = nullptr;
  float *db = nullptr;
  float *dc = nullptr;
  cudaMalloc((void **)&da, size);
  cudaMalloc((void **)&db, size);
  cudaMalloc((void **)&dc, size);
  cudaMemcpy(da, ha, size, cudaMemcpyHostToDevice);
  cudaMemcpy(db, hb, size, cudaMemcpyHostToDevice);
  cudaMemcpy(dc, hc, size, cudaMemcpyHostToDevice);
  
  // kernel config
  const int thread_per_block = 256;
  const int block_per_grid = (size + thread_per_block - 1) / thread_per_block;
 
  VectorAddGPU<<<block_per_grid, thread_per_block>>>(da, db, dc, size);
  
  // 把数据从GPU拷贝回CPU
  cudaMemcpy(hc, dc, size, cudaMemcpyDeviceToHost);

  // 释放GPU显存
  cudaFree(da);
  cudaFree(db);
  cudaFree(dc);
    
  // 释放CPU内存
  delete[] ha;
  delete[] hb;
  delete[] hc;
    
  return 0;
}

deformable conv
- torch已经有官方版本了：https://pytorch.org/vision/stable/_modules/torchvision/ops/deform_conv.html
- 民间版本：https://github.com/chengdazhi/Deformable-Convolution-V2-PyTorch/tree/pytorch_1.0.0
- 也是weighted bilinear这部分需要cuda kernel实现
- src下面是cuda files
- functions下面是pytorch的function封装，相当于torch.nn.functional.conv2d这样使用
- modules下面是pytorch的module封装，相当于torch.nn.Conv2d这样使用
MSDeformAttn
- kernel forward & backward在ops/src/cuda/ms_deform_attn_cuda.cu中
- 打包在ops/src/vision.cpp中
- 编译在setup.py中
- torch op封装在ops/functions中MSDeformAttnFunction，用静态方法定义pytorch op的forward & backward

python bindings

从 Python 调用 C 或 C++，以及将数据从 Python 传递到 C 或 C++
ctypes直接加载dll/so

PyBind11

PyBind11 生成的 Python 绑定是一个完整的 Python 模块，可以直接导入和使用

step1：在c++侧编写方法

// cmult.h
float cmult(int int_param, float float_param);

// cmult.cpp
float cmult(int int_param, float float_param) {
    float return_value = int_param * float_param;
    printf("    In cmult : int: %d float %.1f returning  %.1f\n", int_param,
            float_param, return_value);
    return return_value;
}

step2：pybind module编写绑定
1
2
3
4
5
6
// pybind11_wrapper.cpp
#include <cmult.h>

PYBIND11_MODULE(cmult_example, m) {
m.def("cmult", &cmult, "A function that multiplies two numbers");
}
PYBIND11_MODULE宏是python解释器导入扩展模块的入口，当调用自定义的python库函数cmult时，cmult的函数地址的将被调用

step3：setup编译并添加python module

# setup.py
from pathlib import Path
from pybind11.setup_helpers import Pybind11Extension, build_ext
from setuptools import setup


example_module = Pybind11Extension(
    'cmult_example',
    [str(fname) for fname in Path('src').glob('*.cpp')],
    include_dirs=['./'],
    extra_compile_args=['-O3']
)

setup(
    name='cmult_example',
    description='pybind11+setup.py example',
    ext_modules=[example_module],
    cmdclass={"build_ext": build_ext},
)

安装：python setup install，就会添加名为cmult_example的python module

step4：python脚本中调用module以及cmult函数

from cmult_example import cmult

a = 10
b = 1.25
print(cmult(a, b))

Cython

pytorch-quantization

tensor quant和fake tensor quant：given a tensor，调用pytorch_quantization.tensor_quant.tensor_quant(x)返回一个tensor的量化tensor(int)和scale，调用pytorch_quantization.tensor_quant.fake_tensor_quant(x)返回一个tensor的伪量化值(float)，即QDQ的值，可以对比原始tensor，计算quant error
QuantDescriptor：量化配置
TensorQuantizer：基于量化配置实例化一个Quantizer

quant_nn：实例化一个quant layer

given fc1 = nn.Linear(in_features, out_features, bias=True)

quant_fc1 = quant_nn.Linear(
    in_features, out_features, bias=True,
    quant_desc_input=tensor_quant.QUANT_DESC_8BIT_PER_TENSOR,
    quant_desc_weight=tensor_quant.QUANT_DESC_8BIT_LINEAR_WEIGHT_PER_ROW)

c++ basics recall
- cmake
- 编译原理
- 堆栈
- L1/L2 cache
- 深浅拷贝、传变量/传引用
- 重载：函数名相同，用参数和功能来自动区分
- 类和对象
  - 类：抽象，定义属性，默认private，protected只能被子类访问，public可以被访问
  - 对象：类的实例，具有类的所有变量和属性
- 继承
  - class derive_cls：public parent_cls1, public parent_cls2
  - 多态：重新实现父类的方法，来实现不同功能
自动驾驶 basics
- bev：BEVFormer
- 3D det：PETR
- 多模态：lidar+camera

# raw conv forwards
def conv2d(X, W, B, stride, pad, dilation=1):
    # X: [b,c,h,w]
    # W: [k1,k2,in,out]
    # B: [out]

    # pad
    pad_h, pad_w = pad
    X_pad = np.pad(X, [[0,0],[0,0],[0,pad_h],[0,pad_w]])
    b, _, h_pad, w_pad = X_pad.shape

    # out shape
    s, d = stride, dilation
    k1, k2, _, out_channels = W.shape
    k1, k2 = k1+(k1-1)*(d-1), k1+(k2-1)*(d-1)     # receptive field
    out_rows, out_cols = (h_pad + pad_h - k1) // s + 1, (w_pad + pad_w - k2) // s + 1

    # tranverse spatial
    Y = np.zeros((b,out_channels,out_rows, out_cols))
    for i in range(out_rows):
        for j in range(out_cols):
            i0,i1 = i*s, i*s + k1
            j0,j1 = j*s, j*s + k2
            Y[:,:,i,j] = X_pad[:,:,i0:i1:d,j0:j1:d].dot(W.T) + B  # [b,out]

    return Y


# raw conv backwards
def conv2d_bp(X_pad, W, B, Y_grads, stride, pad, dilation=1):
    # X_pad: [b,c,h_pad,w_pad]
    # W: [k1,k2,in,out]
    # Y_grads: [b,out,h_out,w_out]
    # return: X_grads, W_grads, B_grads

    b, _, out_rows, out_cols = Y_grads.shape
    s, d = stride, dilation
    k1, k2, _, out_channels = W.shape
    k1, k2 = k1+(k1-1)*(d-1), k1+(k2-1)*(d-1)

    # compute grads
    X_grads = np.zeros_like(X_pad)
    W_grads = np.zeros_like(W)
    B_grads = np.zeros_like(B)
    for i in range(out_rows):
        for j in range(out_cols):
            i0,i1 = i*s, i*s + k1
            j0,j1 = j*s, j*s + k2
            kernel = Y_grads[:,:,i,j]   # dL/dY
            B_grads += kernel * 1       # dL/dY * dY/dB
            W_grads += kernel * X_pad[:,:,i0:i1:d, j0:j1:d]   # dL/dY * X
            X_grads[:,:,i0:i1:d, j0:j1:d] += kernel * W       # dL/dY * W

    # remove pad
    pad_h, pad_w = pad
    b, _, h_pad, w_pad = X_pad.shape
    X_grads = X_grads[:,:,:h_pad-pad_h, :w_pad-pad_w]

    return X_grads, W_grads, B_grads

im2col & col2im
- conv2d里面X_pad[:,:,i0:i1:d,j0:j1:d].dot(W.T)这个做的是patch[b,in,k,k]和filter[out,in,k,k]的点积
- im2col可以将循环转化成矩阵乘法
  - filter转化成[out, in*k*k]，X转化成overlap patches[in*k*k, b*n_rows*n_cols]
  - 结果得到[out, b*n_rows*n_cols]
- col2im转换矩阵结果
  - [b, out, n_rows, n_cols]
BN/LN/GN/IN
sgd
adamw
regularization

bit hacks：https://zhuanlan.zhihu.com/p/37014715

判断奇偶：与1做交
1
2
def isOdd(n):
return n & 1

交换两个数字：

def swap1(a,b):
    a = a + b       # a+b
    b = a - b       # a
    a = a - b       # b

def swap2(a,b):
    a = a^b
    b = a^b         # b=(a^b)^b=a^(b^b)=a^0=a
    a = a^b         # a=(a^b)^a=b^(a^a)=b^0=b

将比特序列反向：一个右移弹出，一个左移加入

def inverse(n):
    result = 0
    while n:
        tmp = n & 1
        n >>= 1
        result += tmp
        result <<= 1
    return result

对2的幂数取余

1
2
3

# d=1<<s, {1,2,4,...}
def residual(n, d):
	return n&(d-1)

计算log2：寻找最高位

def log2(n):
    cnt = 0
    while n>>1:
        n >>= 1
        cnt += 1
    return cnt

计算log10：换底数，log10(x) = log2(x) / log2(10)， log2(10)=1233/4096

def log10(n):
    PowersOf10 = [1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, 1000000000]
	t = (log2(n)+1) * 1233 >> 12
    r = t - (v < PowersOf10[t])

判断正负（零）

1 2	sign = -(v < 0) sign = (v > 0) - (v < 0)

判断两数异号

1	sign = ((x ^ y) < 0) // 有且仅有一个负数flag

判断是否是2的幂数
1
sign = n & (n & n(n-1))==0

计算bit数中1的个数

# 计算效率与n的bit位数成正比
def cnt_1(n):
	cnt = 0
    while n:
        tmp = n&1
        cnt += tmp
        n >>= 1
    return cnt

# 计算效率与n中1的个数成正比
def cnt_1(n):
	cnt = 0
    while n:
        cnt += 1
        n &= n-1       # remove the highest 1
	return cnt

计算平方根