CUDA 学习分享

Last Updated on 2023年3月28日

Some small talk

总的来说,使用CUDA和写汇编是比较相似的: 你总是需要知道你写的Cpp代码大致会被编译器生成出什么样的PTX/SASS, 同时大概了解机器是如何执行这些指令的(Have a accurate mental model you are using), 才有可能写出高性能的代码

编程模型核心的三部分

  1. Thread-hiearchy: 从Grid到Thread, 中间存在很多逻辑层,这些逻辑层是如何对应到硬件上的,如何被硬件执行。
  2. Memory-hiearchy: 硬件上可以使用的存储单元有哪些,这些单元的基本特性,容量、带宽、约束,如何使用这些单元。
  3. 同步控制: 不同层级如何进行同步,同步的开销如何。

这三部分都是在学习过程中应当不断关注的. 这篇文章将分Misc和Optimization两部分,Misc主要是一些零碎的知识点,Optimization则是一些优化方法的经验总结.

本文将会假定你已经基本阅读过CUDA Programing Guide, 因此不会在基础的内容上过多着墨.

一个比较合适初学者的学习路线是 "Programing Guide" -> "Professional CUDA C Programing Guide" -> "Tuning Guide" -> "Cuda Handbook 2nd" -> "XXX Whitepaper"等

  • 这样一轮学习后,你就已经基本入门CUDA了,最起码可以听懂别人在说什么.经过一些实践/代码阅读后,在大部分场景下都能写出性能不错的代码
  • 尽管"Professional CUDA C"已经是Kepler时代的资料, 书中的一些基础理念仍然是非常有价值的。

Miscs

Compilation

  • Installation, 参考 https://developer.nvidia.com/cuda-downloads 即可. 唯一的建议是仅安装toolkit, 而不要安装驱动, 以避免驱动更新导致新的问题.
    • 使用apt安装时,不要执行apt install cuda,而是apt install cuda-toolkit-X.Y
    • 使用runfile安装, 在弹出的提示中, 可以选择不安装驱动
    • 如果需要安装驱动, 建议使用已经发布了至少半年的kernel.
    • 另外, CUDA Toolkit 安装完成后, 默认不会添加到PATH中, 相关的环境变量可能仍然需要手动设置.
  • Core types of Driver API:
    • 这些type都是handle-based opaque type
    • cuDevice: 概念上对应一个物理GPU
    • cuContext: 概念上相当于进程, 一个Device可以有多个Context
    • 每一个host thread都有一个 context-stack, 栈顶的context作为implicit used context, ,只有当栈顶对应的context处于激活状态时,thread内的Driver API调用才是有效的.
    • 每一个device都只能有一个激活context, context switch的开销很大.
    • context在host逻辑上是一个shared_ptr,并不是unique_ptr,这意味着context可以在线程间共享, push,attach等操作将会使引用计数增大, detach,destroy等则会使引用计数减小.
    • cuModule: 概念上相当于DSO, 一个Context可以动态加载多个Module, 和DSO类似, 加载过程中,也会涉及DSO的ctor和dtor,例如, ctor阶段会为DSO使用的global mem分配空间,并设定相关的初始值.
    • 拥有全局语义的对象，如__device__, 逻辑上生命周期和Module一致, 会在(Module)被加载时创建,在Module销毁时被销毁, 且Module每加载到不同的Context,都会有不同的实例.
    • cuFunction: 概念上就是Function, 可以用于产生kernel-launch, 可以用类似dlsym的方法从当前context中查找
    • cuEvent: 和cudart的Event一致
    • cuStream: 和cudart的Stream一致
    • cuDevptr: 代表指向device地址的指针
  • 常用的C++ level mark:
    • __attribute__((const)): 函数执行过程不会修改任何外部状态(no side effect), 计算结果也不依赖任何global (return value only depend on argument)
    • __attribute__((pure)): 函数执行过程不会修改任何外部状态(no side effect), 计算结果可以依赖global (return value depend on argument and global states)
    • __CUDA_ARCH__: 只应该用在device相关的code
    • maxrregcount 和launchbound是控制寄存器用量的主要手段, 进而影响occupancy
    • __grid_constant__ 用于修饰kernel的参数, 这种参数会始终保留在const memory中, 而不会拷贝到每个thread的stack上(local-memory),以避免额外开销,也因此对每个thread而言,都有同样的地址.
    • 启用__restrict__可能导致大量的内存访问都被mem2reg优化为寄存器访问, 进而导致occupancy降低
    • 另外注意, 要避免alias的话，所有pointer参数都要有restrict标记
    • builtin_assume, builtin_assume_aligned, builti_exxpect,builtin_unreachable 可以辅导编译器生成代码
  • nvcc编译的所有文件都会在头部自动#include一系列头文件,如<cuda_runtime.h>,这可能导致一些超出你预期的行为.
  • CUDA中有virtual architecture的概念，用compute_xy描述，用于描述PTX版本； 还有real architecture的概念，用sm_xy描述， 对应硬件架构， 实际的编译过程中， 总是会先生成PTX汇编，再基于PTX汇编生成SASS， 有很多编译参数可以控制PTX和SASS使用的版本，实践中建议使用–gen-code,它的歧义最少,且能使用多次

下面的编译选项中,最终的fatbin中将包含3个cubin文件和一个ptx文件

nvcc x.cu
-gencode arch=compute_50,code=sm_50 # 按compute_50指令集输出PTX,然后针对SM_50生成bin,最终只保留bin
-gencode arch=compute_60,code=sm_60 # 按compute_60指令集输出PTX,然后针对SM_60生成bin,最终只保留bin
-gencode arch=compute_70,code=\"compute_70,sm_70\" # 按compute_70指令集输出PTX,然后针对SM_70生成bin,最终同时保留PTX和BIN

注意,这样生成代码时,需要使用编译器生成的host代码才能自动在运行期选择合适的cubin,手动使用driver API只能load特定的cubin,不能直接load fatbin。如果要使用driver API, 应当直接生成单独的cubin文件,然后再用driverAPI去load

  • nvcc自身可以是64位的,也可以是32位的。
    • 当64位时,默认生成64位的device code,当是32位时,默认生成32位的device code。
    • 使用-m64或-m32可以控制device code 使用的默认字长。
    • device code和 host code 的位数必须一致, 否则launch阶段host的数据不能很好的map到device上 i.e., 64位的device code和32位的host code不兼容
  • 在host进程的init,dtor阶段,cuda自身会自动执行一些代码,因此,不应该在init、dtor阶段插入用户的cuda代码,避免UB
  • 由于一定的实现问题, cudart最好是以静态链接的形式链接进入程序, 这一方面便于分发,另一方面可以避免复杂的兼容性问题.
    • CUDA 6 之后,所有CUDA Toolkit自带的库默认都是静态链接的,且总是会静态链接cudart
  • PTX: 总的来说, PTX和nvvm和cuda C基本是平级的,只是多了一些更精细的控制,仍然不能对应到最终的指令.
    • PTX 是 load store machine, 且仍然使用虚拟寄存器
    • CUDA C 的所有概念都在PTX里有直接对应
    • PTX暴露出了更精细的内存控制, 所有的IO指令都可以指定内存空间, 以及是否使用cache
    • PTX提供了更多的builtin寄存器, 如%warpId等, 以及更多的内置函数(部分内置函数可能会因CUDA C更新延迟而没有暴露为C intrinsic function).

Execution

  • Block在分配给SM后,就只会由对应的SM执行, warp在分给某个schedular后,只会被这个scheduler调度.
    • Block分配的策略是比较灵活的,会按空闲资源调度. warp的分配则就是简单的round-robin平均分配.
    • scheduler能使用的SM资源也是被平均分割的,因此一个bolck内thread的数量一般至少需要是scheduler_num * warp_size
  • 自Pascal之后,group和cooperative op的概念将变得非常重要. cooperative op 都会涉及一个明显的group对象, 意思是:由这个group内的所有线程合作完成某个功能(而不是各个线程独立的各做一样的操作).
  • cooperative op 物理上是有两种实现模式:
    • 每个thread会各自负责大任务内的一部分子任务。
    • 由某个thread执行, 其余thread都是no-op
  • 每一个device都有一个唯一的NULL Stream, 它是默认的stream模式(legacy), NULL stream的行为很复杂, 建议一定要使用手动的stream
  • 隐式同步:若调用foo()会触发隐式同步,则在同一个Context中, foo()之前的所有stream都必须执行完成后,foo()对应的操作才能开始执行.
    • 任何发送到NULL Stream的API都会造成隐式同步
    • 注意, 很多没有指明Stream的API都会发往Default Stream, 如果DefaultStream是Null Stream的话,就会造成隐式同步.如cudaMalloc,cudaMallocHost,cudaMemset,cudaMemcpy
    • 各类系统配置类的API也可能触发隐式同步,例如L1相关的参数配置.
  • Programmatic Dependent Launch and Synchronization. Hooper之后, 单个Stream内的相邻block可能交错执行, 例如,两个block有顺序依赖,前一个block在执行write阶段 ,后一个block的prepare阶段就已经可以开始了, 这需要显式的编程实现.
  • "Warp Occupancy": 用于表征一个kernel 的warp level 可并行能力的指标,它的值为(kernel在一个SM内至多的warp数) / hw_max_reside_warp
    • hw_max_reside_warp=hw_max_thread_num/hw_warp_size
    • 例如,一个kernel在launch时产生了8个Block,每个Block内有4个warp. 假如每个block都需要48KB的SMEM,且SM仅有96KB的SMEM, 则一个SM内至多也就能放2个这样的block, 也就对应至多8个warp. 假设这个SM至多支持32个reside warp,那么这个kernl launch也就只能达到 8/32 = 25% 的Occupancy.
    • 类似的概念还有SM Occupancy, 需要通过足够的block来保证
  • 一种测试occunpancy是否充足的方法是:
    • launch时额外动态分配shared mem,使得occupancy降低.
    • 如果降低occupancy后性能没有下降,一般也就意味着occupancy已经足够了.此时,提升occupancy不会带来明显的收益.
  • GPC 对应了 TBC , SM 对应了 Block, TBC内的Block保证分配到同一个物理GPC中, 同一个Block内的Thread保证分配到同一个SM中.
    • TBC的主要特性是支持Distributed Shared Memory, 可以直接获得其他Block内smem的地址
    • 目前只有read/write/atomic op 可以直接用在DSMEM中
  • CUDA Event逻辑上是一个在stream内的一个Op， 当Stream执行到Event时，它就会被set.
    • 每次调用Record时，都相当将Event移动到新stream当时的末尾(也就是说，先前调用的Record就无效了)
    • 当Event被触发后，event内会存储触发的时间
    • Event可以作为一个节点，用于同步多个Stream，例如，下面就形成了一个菱形节点，上下分别是AD，中间是BC

从实现上说, Event的Recrod实际是对当前的Stream做一个快照,这个快照信息将用于决定Event是否完成.

kernel_A<<<  ...,  stream1  >>>(...);
// Record event1 as the "Kernel A Finished"
cudaEventRecord(event1,  stream1);

// just run B after A in stream1
kernel_B<<<  ...,  stream1  >>>(...);

// let stream2 wait kernel_A
cudaStreamWaitEvent(stream2,  event1);
// launch C in stream2, C will start after event1
kernel_C<<<  ...,  stream2  >>>(...);
// Record event2 as the "Kernel C finished"
cudaEventRecord(event2,  stream2);

// Join stream2 back to origin stream (stream1)
cudaStreamWaitEvent(stream1,  event2);
kernel_D<<<  ...,  stream1  >>>(...);
  • cudaStreamCreateWithPriority(&st_high, cudaStreamNonBlocking, priority_high); 可以调整Stream的优先级
  • 逻辑上的并行仍然可能受到物理设备的限制,必要时,应当清楚的知道物理上并行是如何达成的. 例如，逻辑上不同Stream的grid是并行的，但是物理上可能会把它们对应的grid放到同一个物理的Connection上，从而降低了并行性.

Warp Level Instruction and Independent Thread Scheduling:

  • Warp Level Instruction是指由warp内的线程协作完成的指令.
  • Volta之前, 由于隐式同步, warp内的32个线程都必须convergent的参与同一个warp操作
    • 换言之, Volta之前的xxx(args...)等价于Volta之后的xxx_sync(0xffffffff,args...)
  • Volta之后, warp内的32个线程可以只有部分线程参与, 通过mask标记即可, 例如__syncwarp(0x0000ffff)可以仅同步warp内的线程0~15
  • xxx_sync(mask,args...)的意义为: 只有当mask对应的thread都到达同步点之后,才能开始执行xxx操作.
    • 注意: mask仅描述参与同步的线程, xxx操作理论上可以涉及未参与同步的线程(如果这样做,结果显然将是UB). 所以mask称为sync_mask更合适,它对xxx自身的语义没有影响.
  • 从实现逻辑上说, 每个warp都有一个共享的warp_sync_mask, 用于标记到达同步点的线程,当线程执行__xxx_sync(expect_mask)时,效果相当于 
    set(warp_sync_mask , thread_id % 32) // set bit to 1
while(warp_sync_mask != expect_mask){
sleep()
}
clear(warp_sync_mask,thread_id % 32) // clear bit back to 0
__xxx()
  • 显然,只有当expect_mask中的所有线程都到达同步点的时候,才可能放行所有线程. 理解了这个实现, 就能明白CUDA官方文档中对expect_mask的约束
    • 约束1: 所有参与同步的thread都必须体现在mask中.
    • 约束2: 所有参与同步的thread都必须使用同一个mask.
    • 约束3: 所有参与同步的thread都必须使用同一种指令.
  • 这种实现逻辑还带来了一个新福利: 不同branch的同种指令也可以进行同步, 对于下面的例子:
    • 假设前15个thread先执行,那么它将到达shf1并hang住,调度器切换到后15个thread开始执行,后15个thread将到达shlf2,此时,mask对应的30个thread都到达了同步点, 同步完成后,才开始执行shfl操作.
    • 另外,这里, shfl1和shfl2不但可以进行同步, 各自的参数也是不同的
    • 相比之下,Pascal的实现就显得比较复杂了.
#include <cassert>
#include <iostream>
constexpr int sync_thread_num = 30;

// Only valid after volta
__global__ void strange_shuffle_volta() {
  assert(blockDim.x == 32);
  int swapped = 188;

  // only thread 0 to 29 join the shfl
  constexpr uint32_t mask = 0x3fffffff;

  if (threadIdx.x < 15) {
    int val = threadIdx.x;
    swapped = __shfl_sync(mask, val, threadIdx.x + 1); // shfl1
  } else if (threadIdx.x < 30) {
    int val = threadIdx.x * 2;
    swapped =
        __shfl_sync(mask, val, (threadIdx.x + 2) % sync_thread_num); // shfl2
  }
  printf("thread %d get %d\n", threadIdx.x, swapped);
}

// One possible implementation before volta
__global__ void strange_shuffle_pascal() {
  assert(blockDim.x == 32);
  int val = 0;
  if (threadIdx.x < 15) {
    val = threadIdx.x;
  } else if (threadIdx.x < 30) {
    val = threadIdx.x * 2;
  }
  // Note: __syncwarp is not necessary for pascal ,this is just to make sure
  // running correctly on later SM
  __syncwarp();
  int swapped_low = __shfl(val, threadIdx.x + 1);
  __syncwarp();
  int swapped_high = __shfl(val, (threadIdx.x + 2) % sync_thread_num);

  int swapped = 188;
  if (threadIdx.x < 15) {
    swapped = swapped_low;
  } else if (threadIdx.x < 30) {
    swapped = swapped_high;
  }
  printf("thread %d get %d\n", threadIdx.x, swapped);
}

int main() {
  printf("strange_shuffle_volta\n");
  strange_shuffle_volta<<<1, 32>>>();
  cudaDeviceSynchronize();
  printf("strange_shuffle_pascal\n");
  strange_shuffle_pascal<<<1, 32>>>();
  cudaDeviceSynchronize();
  return 0;
}
  • 最后,__activemask()返回的mask是另一回事, 它仅仅标明当前的激活线程. 从理论上说,它是不可预测的, 所以几乎没有实用价值. 例如,对于下面的代码,编程模型允许warp Scehduler分32批执行warp内的线程, 也就是说, 先执行第一个线程,打印出0b1=1,再执行第二个,打印出0b10=2. 
    __global__ void foo(){
print("%d",__activemask());
}

Dynamic Parallelism

将host侧的 runtime API 都porting到device侧, 进而支持在device上launch kernel, 就是dynamic parallelism的主题. device侧的API行为和host一致

  • 逻辑上:

    • grid内的任意thread都是对等的,相当于host侧的thread
    • stream，event等在device thread上创建的对象逻辑上grid private的，只有创建的grid能使用它们.
    • host， grid， subgrid是三个scope，三者之间的runtime对象不能共享，如stream，event

  • device侧的impilict stream，null steram行为也很奇葩，最好不要使用

  • 在laucnh sub grid的时, 只保证subgrid见到的内存和当时的thread是一致的, 这意味着一般需要在launch之前加入一些level更高的memfence，或者sync，以保证launch时能看到一致的内存

  • Stream

    • The Fire-and-Forget Stream 字面意思： 就是创建一个临时stream，再issue进去， 这个语法糖硬件上更快， The Fire-and-Forget Stream 总是独立的，不能被任何手段sync.
    • TailStream: 保证新kernel在当前grid结束后才开始执行
    • 假如某个stream初始只有[grid0,grid1],那么实际的执行顺序是: {grid0 and grid0’s none tail sub-grid} {grid0’s tail graid} {grid1 and grid1’s none tail grid} {grid1’s tailgrad}
    • gridx and gridx’s none tail sub-grid 彼此都是并发的
    • gridx’s tail graid 是串行执行的,按插入顺序launch

  • cuda graph也可以在device侧launch. 但必须issue到 cudaStreamGraphFireAndForget，cudaStreamGraphTailLaunch之一的stream，

    • 在device上调用cudaGraphLaunch来launch graph的kernel, 其Dynamic Pallasim是被禁用的, 因此不能 launch sub-grid，

  • 如果某个kernel会launch graph, 这个kernel也可以作为node添加到另一个graph中，这将可以形成graph的嵌套.

barrier

barrier的工作逻辑和cpp完全一致

  1. 需要有一个唯一的对象被参与同步的所有线程共享(存放在shared或device中),且这个对象只能被初始化一次,所以类似下面的代码总是必须的
  __shared__  cuda::barrier<cuda::thread_scope_block> bar;
  if (group.thread_rank() == 0) {
      init(&bar, group.size());
      // Single thread initializes the total expected arrival count.
  }
  group.sync()
  1. CUDA对barrier的latchdown做了更多的限制: 调用wait时,token的phase至多只能比当前bar.phase慢1,否则将是UB, 也就是说,下面的代码将是UB
// bar是一个block_size的barrier, block内每个thread都执行
auto token = bar.arrvie(); // phase advanced
bar.arrive();// phase advanced again
bar.wait(std::move(token)) // UB

barrier内部会持有count和phase_id. arrive()是执行latchdown的核心操作,当latchdown到0时,会先执行completionFn,再尝试放行所有wait的线程,最后返回当前的pahse_id,再重置count,更新phase_id的值。

arrive会返回一个token, wait(std::move(token))意思是等待token对应的phase结束

cuda::barrier有一个thread_scope的模板参数,该模板参数主要用于描述共享的范围,没有任何其他意义,

  • thread_scope_system,thread_scope_device,thread_scope_block是三个可选项,默认的cuda::std::barrier等价于cuda::barrier<thread_scope_system>
  • 逻辑上说, thread_scope_device的barrier也可以存储在 shared mem 内, 当做thread_scope_block的barrier用, 不过这样性能会变差.

虽然功能上说barrier可以用于同步所有的线程,不过在这些场景,使用__syncthreads, group.sync()之类会更加高效。

memcpy_async

memcpy_async的API有三个变形, 分别是group,barrier,pipeline.

  • group: 搭配cooperative_groups::wait(group)使用,注意不是 group.sync(), 从逻辑上说，相当于动态为group插入了一个异步线程.cooperative_groups::wait只是让当前线程等待所有异步线程完成,并不会有group.sync()的效果.
  • barrier: 搭配barrier.arrive_and_wait()使用,从实现角度讲,相当于动态增大了barrier的count,这个count只有在async_memcpy完成时才会下降
  • pipeline: 约等于线程安全的queue

注意, memcpy_async的变形很多,使用时一定要注意API对应的group, 因为这直接决定了memcpy时,src/dst/size的值.

pipeline

它是一个专门针对软流水的语法糖, 当需要实现软流水时,总是应该优先使用它(虽然使用barrier也可能实现相似的功能)

    __shared__ cuda::pipeline_shared_state<
        cuda::thread_scope::thread_scope_block,
        stages_count
    > shared_state; // shared_state is core object in shared memory
    auto pipeline = cuda::make_pipeline(block, &shared_state); // pipeline is a proxy object on the thread stack

  1. 它实现上是一个proxy pattern, cuda::pipeline是每个thread访问pipeline_shared_state的proxy

    • pipeline_shared_state的模板参数也仅仅是描述pipeline会被共享的范围,和barrier类似.
    • pipeline_shared_state需要在共享的内存区域创建
    • thread_scope的pipeline是性能最优秀的, 它不使用任何共享资源, 用cuda::pipeline<cuda::thread_scope_thread> pipeline = cuda::make_pipeline()直接创建
    • make_pipeline是一个同步操作,它用于初始化pipeline_shared_state,为当前线程确定role, 并通过线程间通信确定group内producer/consumer的数量

  2. pipeline逻辑上是一个fifo, head in, tail out, 这个pipeline的元素称为stage, 它逻辑上和git的staging area类似,用于暂存memcpy_async等指令.

    • pipeline这个proxy有三种可能的角色, consumer, producer, both
    • fifo的最大容量是编译期创建pipeline_shared_state时指定的,当fifo中的stage满时,后续的producer将在acquire时被阻塞.
    • 对于role为producer的pipeline:
      • pipeline.producer_acquire(); 表明当前thread开始尝试向fifo中push数据,并需要lock相关的资源
      • 当producer_acquire成功后,当前thread的后续async指令都可以发射到acquire获得的stage中
      • pipeline.producer_commit();提交当前thread的async任务
      • 当group内所有的producer都commit后, 对应的stage将正式被push到fifo中,这个stage将会在内部所有async任务完成后被标记为ready
    • 对于role为consumer的pipeline
      • pipeline.consumer_wait(); 表明当前thread需要从fifo中取数据,并需要锁定相关资源
      • 当consumer_wait()成功后,说明group中的producer提交的一个stage已经ready,可以开始处理数据
      • pipeline.consumer_release(); 表明当前thread的任务已经完成.
      • 当group内所有的consumer都release后,这个stage将被pop出fifo
    • 注意: stage只有逻辑意义, fifo tail对应的是哪些数据需要由用户自行track

Memory

  • cuda device malloc 使用的是独立的device分配器,使用的是一段称为device heap的内存段( which may be smaller than the available unused device memory),只能在device侧动态分配/销毁
    • cudaMalloc() and cudaFree() When invoked from the device, these functions map to device-side malloc() and free().
  • 在sm35之后, Cache Line大小基本都是128B, 所以thread的单次IO最好都是以4K byte为单位, 以尽可能保证warp对应的数据对齐到128B
    • 4的倍数越大越好, 例如使用int4等.
    • 可以使用half2这样的凑足4B
  • 对managed memory, Pascal以后的架构支持了Page Fault机制, Host端的主存从逻辑上变为了Device memory的Swap, 是以Page为单位交换的,而不是整段内存,所以效率变高了, 同时,内存的oversubscription也可以被支持
  • CUDA array 主要是为了texture feching优化的(推测是有特殊的fizzle格式)，也就是说，其取值的方法主要是fetch(x_pos,y_pos)，中间会涉及到插值，边界clamp的问题， 而不是按照顺序索引取值
  • set-aside是一个对L2部分控制的机制，在CC8.0以上可以支持
    • 大概率会重复访问的数据被称为Persisting， 只使用一次就再也不使用的数据被称为Streaming， 这两者显然对应了不同的Cache风格。
    • 具体而言,需要
      • 将L2的一部分set-aside为Persisting模式
      • 将Stream内kernel会访问的某一部分内存片段标记为Persisting，意思是在这个Stream内的Kernel执行时，对应的内存片段会被重复访问，例如，前一个Kernel的Output被后一个Kernel的input使用时，相关的内存段就可以被标记
    • 标记为Persisting的内存段将优先由Persiting的L2服务
    • L2中的set-aside区域是所有Stream共享的， 也就是说， 如果跨Stream的并行程度很高，那么仍然可能有性能问题。
  • set-aside使用过程中主要由三个参数影响： L2的预留尺寸 site-aside-size, 被标记的内存片段window-size， 以及hit-ratio:
    • 硬件会按固定pattern选择单个内存窗口内 window-size * hit-ratio 的数据暂存在site-aside-size的缓存中。
    • 这也就意味着，当所有并行Stream的 window-size * hit-ratio 之和大于 site-aside-size 时， 后缓存的数据仍然会evict先前已缓存的数据， 甚至可能产生抖动
  • Pinned Memory:
    1. Portalble: 在没有启用UVA时,pinned Memory只对当前的device有效,需要在mallochost的时候额外传入cudaHostAllocPortable对所有设备有效,启用UVM时,内存自动是portable的
    2. 在mallocHost时,cudaHostAllocWriteCombined, pinned memory 上可以禁用Cache,这将导致host侧读取性能变得非常差,但是相应的,device 侧的IO性能可以变得很好。
    3. Mapped,可以额外传入参数, 把Host memory映射到Device的地址空间,从而可以直接让kernel自动穿透PCI-E访问host侧的数据。在启用UVA时, 内存自动是mapped。
  • kernel-launch时, 参数会作为一整块buf从host拷贝到device的constant memory, 这就从原则上要求kernel的参数必须是trivially copyable and trivially destructible.
    • 还会受到const memory cache 48KB 的限制.
  • Unified memory address未启用时, device对象在device侧的地址和host侧的地址可能是不同的
    • 例如对于一个全局变量__device__ int a;,它的地址可以在device侧通过&a获得,也可以在host侧通过GetSymbolAddress()获得, 二者可能不同.
    • 这就意味着: device侧获得的地址只能在device代码使用, host侧获得的地址只能在host侧使用. 例如GetSymbolAddress()获得的地址可以直接用在cudaMemcpy()中,但是不能把这个地址值传入到device侧解引用.
  • 硬件无法保证aliased virtual adress的正确访问。
    • __global__ void foo(char *A, char *B)中，如果A,B的运行期值不同，但是却通过virtual mapping对应到同一个物理地址，那么这种访问将会有一致性的问题。
  • 从硬件的角度看:
    • SM内的4个Scheduler会共享一个ld/st单元,这意味着共享了到L1 (shared memory),L2,Global的带宽.例如, 假如有4个warp占用了4个scheduler,每个warp都issue了一条从shared mem 读取数据的指令, 这4个指令物理上是串行被执行的.
    • 不同的SM到L2,Global的出口是独立的, 可以并行使用. 官方标称的Global <-> REG 带宽, 实际上是由很多个SM组合起来达成的, 标称虽然很大, 但是每个SM分到的其实仍然比较小.
    • 一般来说,对于单个Schedular, 从L1到register的带宽是64 byte/CLK, 它是一个常量. 但是latency则是变量, latency约等于指令byte数的一半. 例如,假如单个thread lds 8byte数据,那么latency就是4CLK, lds 4byte数据,那么latency就是 2CLK.
    • 例如, 当一个warp issue 一条 lds.4 时,那么相当于发出了32 * 4 = 128 byte 的 request, 这个request 将在 2个CLK 之后完成, 相关的register才会被标记为可用.(我们也容易计算出, 等效带宽就是 128 / 2 == 64)

ThreadFence

总的来说,CUDA的各类threadfence和CPP的threadfence工作逻辑是相同的,但是由于cuda的thread_fence()无法影响到L1 cache, 无法保证跨SM的L1一致性,只能在L2, Global三个位置建立正确的全局happens-before关系。

具体而言，编译器针对threadfence仅保证物理上的write/load 的执行顺序，但是并不保证能否取到正确的值。

换言之,CUDA threadfence的visible、visibility是有约束的,并不是无条件成立：

  1. 在block内,happens before 总是成立的(因为不会跨越SM)
  2. 在可能跨越SM时,happens before 需要volatile配合,强制编译器生成绕过L1的指令,才能保证happens before成立

CUDA内所有的同步指令都会带有特定scope内的memory fence。 例如,__syncthreads()会带有一个block level的memfence,对于block内的thread, L1,L2,global 总能保证一致,所以官方文档也会说保证__syncthreads()后,之前的内存都是visibile

Tensor

  • 常见Layout
    • NC/kHWk表示按照两级Tile来遍历，
      • 例如 (1,64,5,4)的数据如果按照NC/32HW32的结构排布，那么相当于tile的阵列为(1,2,5,4)，内层每个tile是(1,32,1,1),然后外层按照NCHW的顺序遍历，内层则是逐元素存储
    • 对于矩阵[B,M,N]，就是RowMajor和ColMajor，对应的Stride分别是
      • Row: [MN,N,1] , 也就是说 A[c,i,j]位于 c*MN+i*N+j
      • Col: [MN,1,M], 也就是说 A[c,i,j]位于 c*MN+i+j*M
    • Packed: 内存完全连续，stride等价于prod
    • Partialy packed: 高dim上可能不连续，低dim上是prod
  • interleave: CUDA支持特定pattern的interleave, 简单来说,就是连续的X个128 byte数据中,以16byte为单元交换位置,每 128 byte有一个独立的子pattern, 每 128X byte为一个周期

Tensor 的im2col 访问模式注释:

ptx_im2col_layout

  • 在im2col中,pixel这个词对应的是一个向量,存储了chanel方向的所有值. 也就是说 (N,D,H,W,C)的图像中,只有 (N,D,H,W)个像素,每个像素有C个channel
  • box
    • corner,和N无关, 是直接划定边界, lower(d0,h0,w0) ,higher(d1,h1,w1), lower是相对左上角的坐标,higher是相对右下角的坐标, 就是 [d0,h0,w0] -> [D+d1,H+h1,W+w1] 对应的box
  • 起点:
    • coordinates, 第一个元素的坐标,注意,是元素,不是像素
    • offset, 用于自动在H和W方向上增加偏移, 这个偏移会影响coordinates和corner
  • 遍历参数:
    • Pixels-per-Column: P,从起点开始, 在box内, 连续遍历P个pixel (把DHW三个方向拉直,也就是说,如果某个W方向访问完了,进入下一个height, HW访问完了,会自动进入下一个depth)
    • channels-per-pixel:C, 遍历pixel的过程中,仅保留每个pixel中从coordinates[-1]开始的连续C个channel的值
    • stride: 在遍历P个像素时,相邻两个像素的步长差

Optimization

  • 名词
    • Inst Thoughput: Scheduler向某个单元 issue 一个指令后, 下一次向同一单元issue的最小间隔. 一般用 issue 1 instruction every K clk描述, 或称为CPI (Cycles per Instruction)
    • Inst Latency: 单条指令的结果变为可用状态(可供后续依赖它的指令使用)的clk数称为Latency，即 complete 1 instuction in L CLK
    • TLP: A schedular issue multiple instruction from different warp
    • ILP: A schedular issue multiple instruction from one warp
  • 指令的 CPI <= Latency.
    • 对于没有流水的单元, CPI 等于 Latency, 也就是说, 每完成一条指令, 才能issue下一条, 内存单元可以认为是完全没有流水的单元.
    • 大部分计算单元都有流水,即便是tensor core,也有compute-write两级流水.
  • 大部分时候, 为了便于分析, 我们都是取一个折中值, 来近似为没有硬件流水的情况. 例如, 对于CPI为1, latency 为6 的指令, 如果程序中大部分计算之间都是没有数据依赖的, 那么我们可以近似用1来作为latency, 也就是每过1个clk,能issue一个新的指令,且上一条指令计算完成.

Tensor Core 中, 可以认为一个clk完成 m*n*k 的矩阵计算. 例如, 对于Volta, 单个clk可以处理的fp16矩阵乘为4x4x4, 那么对于 16x8x8 的矩阵乘, 它可以拆分成 16个 4x4x4, 共需要16个clk. 如果通过单个warp来处理这个16x8x8的矩阵乘, 那么由于Warp有2个Tensor core, 所以需要8个clk才能完成计算, 正因如此, Volta中的 hmma.1688矩阵乘CPI为8. 由于结果还写回寄存器也需要时间, 所以实际的latency会偏大一些, 为14个clk.

注意, 本文说的clk都是schedular的CLK, 因为内存单元一般会有自己的时钟,内存对应的clk一般可以通过延迟时间和schedular的频率计算出来,例如,假如已知内存的latency是89ns, schedular的时钟是1.5GHz, 那么对schedular而言, 内存的latency就是133.33clk

Always recheck Compute Capabilities

ILP and TLP

Given that I[i][j], which means the i’th warp’s j’th instruction. Assume all instructions are fully pipelined, i.e., if there is no dependency, every clk could issue a instruction.

  • e.g 1, TLP
    • Scheduler issued I[0][0] first, for I[0][1] has dependence on I[0][0], I[0][1] can only be issued L clk later, then, scheduler select I[1][0] to issue, and samething happened again, then scheduler switched to warp 2,and so on. To make scheduler busy enough, this schedular should have at least L warp to issue, so that L clk later, the I[0][1] will be ready to issue, and we can restart from warp0 again. If a SM have 4 schedular, it should have at least 4L warp to run, which means the block should have at least 4L*32 reside-threads.
  • e.g 2, TLP and ILP
    • this time, Inst[0][0] and Inst[0][1] has no dependency, so schedular can issue I[0][1] after I[0][0], but I[0][2] depend on I[0][0], so scheduler select another warp again, and this time, schedular only need slect L/2 warps to wait I[0][0] to be completed. And, if a SM have 4 schedular, a block will only need 2L*32 reside-threads this time.
  • In summary: TLP and ILP are equally important technically. Usually ILP is more easy to meet, so ILP is more important than TLP.

Note, 4L ,2L here has no practical meaning, for not all instructions are fully pipelined.

Gneral guidence

  1. Maximize parallel execution to achieve maximum utilization;
    • App Level: consume all device connection by async launch, try not let grids wait each other
    • Block level : consume all SM, and try not let blocks sync each other
    • Warp level: consume all Scheduler， have engough warp, and each warp has enough instruction to issue
  2. Optimize memory usage to achieve maximum memory throughput;
    • Have enough active threads to send memory request in burst style.
    • Minimize transaction in every level
    • Optimize Mem access pattern to get better efficency, choose correct minimal size for each thd, padding if necessary, avoid banck conflict
  3. Optimize instruction usage to achieve maximum instruction throughput;
    • Minimize Low throughput instruction, trade speed to lower precison fp
    • Mniimize divergence
    • Minimize sync
    • Minimize instruction itself
    • Know some cost inst: sp function, int mod, int div
    • Know some fast: half2, not half; rintf(), not roundf(); rint(), not round().

定量分析

Basic Principle

开始定量分析前,需要补充一些基本原则.

  1. Nobdy cares about flops
    • 对于一款合格的硬件, 它的flops总是远远大于带宽的, 问题的难点不是硬件该有多大的flops, 而是软件如何能利用满这么多flops.
  2. Memory is the key, more specifically, the bandwidth and latency of memory is the key
    • 内存带宽是硬件的限制 ,对于内存而言, 最优的访问模式是 burst sequential 模式, 我们需要尽可能发出大的连续内存request, 才有可能达到最大带宽,具体而言,对于带宽为B byte/clk的硬件,假设指令的延迟为L clk,那么应该每L个clk issue 一条请求BL大小的 request, 才能保证带宽始终满载.
    • 可以近似的说, 若 request0 在t0 发出, 那么t0+L时刻, request0 对应的的数据才会被标记为ready.
    • 带宽B是一个硬件约束的常量, 但是L则可能根据指令的不同而改变, 例如, lds.4 byte和lds.16 对应的latency都是不同的
    • 一般而言,推荐使用ld.8, 也就是单个thread每次至少要使用8byte, 每个warp至少要issue 256 byte的request.

通常来说，硬件设计上会使得 B * L 对各级存储都是相近的常量, 这样才能保证同样数量的active-warp有能力充分利用各级的带宽. 例如, DDR是带宽很小Latency很大的单元, L1 Cache则带宽很大Latency很小, 两者应基本满足 B_DDR * L_DDR = B_L1 * L_L1

  1. Usually we dont care about store, because usually we dont wait store

    • 如果真的需要考虑store,那么直接把store串行的放到compute环节之后即可, 也就是说,逻辑上仍然是仅有 load / compute 两条timeline

  2. 我们同样可以通过TLP或ILP来使得IO及计算并行起来, 假设我们严格按照Load,Compute的顺序来组织代码, 那么对于单个scheduler, 至少需要2个reside warp, 才能充分利用所有资源. 例如, 先把 warp0 的 load部分issue出去, 等待L个latency后, 同时把warp0的计算指令和warp1的load指令issue出去, 这样就实现了计算单元和IO单元的并行, 才有可能保证总是有IO指令被执行, 进而到达最大带宽.

    • 也就是说, 在不使用软流水时, 我们至少需要每个schedular都有2个reside warp.
    • 使用软流水时,每个schedular只要有1个reside warp即可.
    • 实践中, 我们应当尽量保证一个block至少有 scheduler * 2 个 warp. 即便使用软流水, 这样也能免费吃到硬件TLP的红利. 对于大部分硬件, scheduler的数量都是4, 所以一个block的thread数至少应该是256.

从内存开始

  1. 确定存储边界

输入数据量为M已知时, 如果某一级存储能完整放下M,那么就以这一级为存储边界.

存储边界的意义为: 我们的优化主要是在存储边界上进行的, 对于更高速的部分, 则不必过多考虑. 例如, 假如我们的代码存储边界是DDR,那么我们在优化过程中,就可以先不必考虑L1, L2, 因为这些部分的带宽都远远大于DDR的带宽, 允许存在一定的低效和浪费.

当我们把存储边界的带宽充分利用之后,才需要考虑下一级的存储,将下一级作为新的边界,继续优化.

SW pipeline

从schedular的角度看, 当我们决定使用软流水时. 我们的指令中要么在issue IO指令, 以load下一次计算需要用的数据, 要么在issue计算指令, 处理已经load好的数据.

假定Load环节有a个独立的load指令, Compute环节有b个独立的计算指令, 那么每个周期的长度就是max(a*Tm,b*Tc), 在完全流水时

  1. a*Tm个clk将用于load下一批计算使用的数据
  2. b*Tc个clk将用于处理上一批load获得的数据
  3. 只有上一批load的数据可用
  4. 只要流水正常, 要么是load的timeline连续,要么compute的timeline连续, 我们的目标是尽可能让compute的time line连续.

虽然我们的理想是让compute连续, 即 b * Tc > a * Tm, 但是这对于很多算法来说是很难的, 因为大部分算法都没有足够的计算强度.例如,对于elt wise指令, 上一批的a个load指令得到的数据, 可能只需要a+C个clk就执行完了, 远小于load下一批数据所需的a*Tm个clk.

对于矩阵乘, 由于它的arith intensity 是线性增长的, 可以近似说, 上一个时间步用 n^2*Tm 个 clk load 的数据, 需要这个时间步中n^3 * Tc 个clk才能计算完成,只要保证n足够大, 想要实现n^3 * Tc > n^2 * Tm还是比较简单的.

最后,从原理上说, 从DDR到计算单元的过程中, 每多一层可编程cache,就能多做一层流水,例如,对于CUDA,有[L0,L1,L2,DDR],其中, L0 register 和L1 shared memory都是可编程的, 所以硬件上CUDA支持2个层次的软流水,我们可以实现最长2x2+1 = 5级的流水线(ddr to smem,smem to reg,compute,reg to smem,smem to ddr).

一般来说,只有三级流水时,结构是(DDR->SMEM,Compute,SMEM->DDR), 这种情况下, 计算密集的算子很容易出现Compute环节很长的情况.
因为Compute环节实际是在操作SMEM,SMEM的带宽及延迟仍然离寄存器有一定距离, 所以Compute环节中可能会有很大一部时间在等待数据从SMEM到寄存器, 这一部分的延迟也是可以被hide的.

single-level
例如,上面的例子中, 如果smem to reg 不进行流水的话, 那么就是lcs,或者lcslcslcs的循环, 但是如果smem to reg 进行流水,那么总的compute时间可以进一步缩短.
mutli-level

在使用多级流水时, 通常会选择

  • rgister: 尽可能大的M作为tiling size, 寄存器中准备2个这么大的buffer,不断从active 的smem buffer中读取数据, 一次处理M大小的数据
  • shared mem: KM作为tiling size, smem中准备两个这么大的buffer,不断从DDR读取数据.

Profiling

当我们选择了足够的并行度, 使用了所有的软流水,编写出代码之后, 理论上说现在的代码就应该免费的撞在roofline上了. 但是一般并不会如此.

此时我们的要做的是: 画出IO及Compute的Timeline, 选择两个timeline中连续的那个作为优化目标,缩短其时间,因为它决定了总时间的长短.

例如, 如果timelime中DDR IO是连续的,那么我们应该优先优化DDR IO, 看是否是存在冗余/不对齐的IO. 如果compute是连续的,那么我们应该优先优化compute.

可能会有一些奇怪的空间时间的tradeoff存在, 例如,将IO的时间缩短一点, 让计算的时间增长一点,可能总的时间会缩短.

在这个过程中, 我们应当不断进行profiling, 判断优化空间还有多少.

  1. 若timeline中DDR连续, 且当前的等效带宽已经接近峰值, 代码层面就没有什么优化空间了,只能换新算法才可能行.
  2. 若timeline中compute连续, 且计算性能已经接近峰值, 那么优化完成.