前言
标准的CUDA编程模型大体上可以抽象为“Prepare&Launch => Execute => Sync => ReadBack”这样的基本流程。在理想情况下，只要计算的workload比较适合数据并行，不存在特别不平衡的划分，并且数据粒度不要太小，这样就能很好地发挥GPU多核的运算性能。但如果底层设备的运算性能非常高，每次完成运算的耗时几乎与kernel launch（~20us数量级）处于相当的水平，那么这种频繁的与host进行交互的工作模式就会带来相对较高的开销，在单个工作线程的场景下会无法打满底层硬件。这种runtime overhead，随着异构计算设备算力的加强，会变得愈发明显。
Persistent Thread（下文简称PT）是一种重要的CUDA优化技巧，能够用于大幅度降低GPU的"kernel launch latency"，降低其Host-Device通讯所带来的额外开销。但由于GPU硬件设计的特性所带来的编程模型限制，导致这种技术没有大规模流行起来（详情请参考综述论文[2]）。本文对PT技术的工作原理，应用场景，代码实现，优缺点以及实现该技术所需要支持的最小功能子集进行分析。
此外，CUDA官方提供了"Dynamic Parallelism"机制用来更好地支持"Data-Dependent Parallelism"这样的特殊情况，所解决的问题与PT技术略有重叠，细节可以参考这个简化版slides以及相应的blog。
Persistent Thread
TL; DR：本质上，PT技术的核心思想就是让同一段kernel代码一直运行在GPU设备上，利用UVA[1]或者zero-copy memory来实现以极低的开销与GPU通信，将workload给offload到计算设备上不断轮询计算任务的kernel，相当于把外部计算设备当做一个CPU可见的线程池。相比之下，通过kernel launch的方式来加载程序（正常的Host-Device交互），需要走一遍完整的从上层runtime一直到PCIE driver的流程，这个开销要大得多。
基本应用场景
CPU-GPU Synchronization
Load Balancing / Irregular Parallelism
Producer-Consumer Locality
Global Synchronization
代码实现
下面以一些具体的代码示例来分析Persistent Kernel这一编程模型的具体应用场景。
设备通讯
考虑如下代码：
// test.cu
// nvcc test.cu -o test && ./test
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#define TIME_INC 100000000
#define INC_CNT 10
#define cudaCheckErrors(msg) \
do { \
cudaError_t __err = cudaGetLastError(); \
if (__err != cudaSuccess) { \
fprintf(stderr, "Fatal error: %s (%s at %s:%d)\n", \
msg, cudaGetErrorString(__err), \
__FILE__, __LINE__); \
fprintf(stderr, "*** FAILED - ABORTING\n"); \
exit(1); \
} \
} while (0)
__global__ void test_kernel(volatile int *data) {
unsigned long time;
for (int i = 0; i < INC_CNT; i++) {
atomicAdd((int *)data, 1);
__threadfence_system();
time = clock64();
while ((clock64() - time) < TIME_INC) {};
}
printf("demo kernel exited\n");
}
int main() {
volatile int *d_data, *h_data;
cudaSetDeviceFlags(cudaDeviceMapHost);
cudaCheckErrors("cudaSetDeviceFlags error");
cudaHostAlloc((void **)&h_data, sizeof(int), cudaHostAllocMapped);
cudaCheckErrors("cudaHostAlloc error");
cudaHostGetDevicePointer((int **)&d_data, (int *)h_data, 0);
cudaCheckErrors("cudaHostGetDevicePointer error");
*h_data = 0;
printf("Kernel starting\n");
test_kernel<<<1, 1>>>(d_data);
cudaCheckErrors("kernel failed");
int value = 0;
do {
int old_val = *h_data;
if (old_val > value) {
printf("h_data = %d\n", old_val);
value = old_val;
}
} while (value < (INC_CNT - 1));
cudaDeviceSynchronize();
cudaCheckErrors("kernel failed");
return 0;
}
上述代码的运行输出为：
$ ./test
Kernel starting
h_data = 1
h_data = 2
h_data = 3
h_data = 4
h_data = 5
h_data = 6
h_data = 7
h_data = 8
h_data = 9
progress check finished
再考虑一个稍微复杂点的Host-Device同步，这里我们用一个progress变量来实时对Host通知kernel的执行进度：
__global__ void
matmult(float *a, float *b, float *c,
unsigned int rowA, unsigned int colA,
unsigned int colB, volatile int *progress) {
unsigned int row = threadIdx.x + blockDim.x * blockIdx.x;
unsigned int col = threadIdx.y + blockDim.y * blockIdx.y;
if ((row < rowA) && (col < colB)) {
float temp = 0.0f;
for (unsigned int k = 0; k < colA; k++) {
temp += a[(row * colA) + k] * b[(k * colB) + col];
}
c[(row * colB) + col] = temp;
if (!(threadIdx.x || threadIdx.y)) {
atomicAdd((int *)progress, 1);
__threadfence_system();
}
}
}
int main() {
int value = 0;
volatile int *d_data, *h_data;
cudaSetDeviceFlags(cudaDeviceMapHost);
cudaCheckErrors("cudaSetDeviceFlags error");
cudaHostAlloc((void **)&h_data, sizeof(int), cudaHostAllocMapped);
cudaCheckErrors("cudaHostAlloc error");
cudaHostGetDevicePointer((int **)&d_data, (int *)h_data, 0);
cudaCheckErrors("cudaHostGetDevicePointer error");
float *h_c, *d_a, *d_b, *d_c;
h_c = (float *)malloc(MAT_DIMX * MAT_DIMY * sizeof(float));
if (h_c == NULL) {
printf("malloc fail\n");
return 1;
}
cudaMalloc((void **)&d_a, MAT_DIMX * MAT_DIMY * sizeof(float));
cudaCheckErrors("cudaMalloc a fail");
cudaMalloc((void **)&d_b, MAT_DIMX * MAT_DIMY * sizeof(float));
cudaCheckErrors("cudaMalloc b fail");
cudaMalloc((void **)&d_c, MAT_DIMX * MAT_DIMY * sizeof(float));
cudaCheckErrors("cudaMalloc c fail");
for (int i = 0; i < MAT_DIMX * MAT_DIMY; i++) { h_c[i] = rand() / (float)RAND_MAX; }
cudaMemcpy(d_a, h_c, MAT_DIMX * MAT_DIMY * sizeof(float), cudaMemcpyHostToDevice);
cudaCheckErrors("cudaMemcpy a fail");
cudaMemcpy(d_b, h_c, MAT_DIMX * MAT_DIMY * sizeof(float), cudaMemcpyHostToDevice);
cudaCheckErrors("cudaMemcpy b fail");
cudaEvent_t start, stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);
*h_data = 0;
dim3 block(16, 16);
dim3 grid(((MAT_DIMX + block.x - 1) / block.x), ((MAT_DIMY + block.y - 1) / block.y));
printf("matrix multiply kernel starting\n");
cudaEventRecord(start);
matmult<<<grid,block>>>(d_a, d_b, d_c, MAT_DIMY, MAT_DIMX, MAT_DIMX, d_data);
cudaEventRecord(stop);
unsigned int num_blocks = grid.x * grid.y;
float my_progress = 0.0f;
value = 0;
printf("Progress:\n");
do {
int value1 = *h_data;
float kern_progress = (float)value1 / (float)num_blocks;
if ((kern_progress - my_progress) > 0.1f) {
printf("percent complete = %2.1f\n", (kern_progress * 100));
my_progress = kern_progress;
}
} while (my_progress < 0.9f);
printf("\n");
cudaEventSynchronize(stop);
cudaCheckErrors("event sync fail");
float et;
cudaEventElapsedTime(&et, start, stop);
cudaCheckErrors("event elapsed time fail");
cudaDeviceSynchronize();
cudaCheckErrors("mat mult kernel fail");
printf("matrix multiply finished. elapsed time = %f milliseconds\n", et);
return 0;
}
上述代码的执行输出为：
$ ./test
matrix multiply kernel starting
Progress:
percent complete = 10.0
percent complete = 20.0
percent complete = 30.0
percent complete = 40.0
percent complete = 50.0
percent complete = 60.0
percent complete = 70.0
percent complete = 80.0
percent complete = 90.0
matrix multiply finished. elapsed time = 6264.967285 milliseconds
事实上在新版CUDA上面，采用cudaHostAlloc分配得到的内存指针可以直接传给kernel就好，这个指针可以通过UVM机制被GPU设备解引用，不需要额外调用cudaHostGetDevicePointer来得到Device可见的指针。这里额外多此一举，只是为阐述其早期编程模型概念。
任务分发
Ping-Pong Buffer
// test.cu
// nvcc test.cu -o test && ./test
#include <stdio.h>
#define TOTAL_ITERS 1000
#define DSIZE 65536
#define nTPB 256
#define BUFFER_EMPTY 0
#define BUFFER_FULL 1
#define cudaCheckErrors(msg) \
do { \
cudaError_t __err = cudaGetLastError(); \
if (__err != cudaSuccess) { \
fprintf(stderr, "Fatal error: %s (%s at %s:%d)\n", \
msg, cudaGetErrorString(__err), \
__FILE__, __LINE__); \
fprintf(stderr, "*** FAILED - ABORTING\n"); \
exit(1); \
} \
} while (0)
__device__ void do_compute(int *buf, int idx, int data) {
buf[idx] = data;
}
__global__ void test_kernel(int *buffer1, int *buffer2,
volatile int *buffer1_ready,
volatile int *buffer2_ready,
const int buffersize) {
static volatile int completed_iter = 0;
static volatile int blkcnt1 = 0;
static volatile int blkcnt2 = 0;
// assumption of persistent block-limited kernel launch
int idx = threadIdx.x + blockDim.x * blockIdx.x;
// persistent until TOTAL_ITERS complete
for (int current_iter = 0; current_iter < TOTAL_ITERS; current_iter++) {
int *buf = (current_iter & 1) ? buffer2 : buffer1; // ping pong between buffers
volatile int *bufrdy = (current_iter & 1) ? (buffer2_ready) : (buffer1_ready);
volatile int *blkcnt = (current_iter & 1) ? (&blkcnt2) : (&blkcnt1);
int my_idx = idx;
while (current_iter - completed_iter > 1); // don't overrun buffers on device
while (*bufrdy != BUFFER_EMPTY); // wait for buffer to be consumed
while (my_idx < buffersize) { // perform the "work"
do_compute(buf, my_idx, current_iter);
my_idx += gridDim.x * blockDim.x; // grid-striding loop
}
__syncthreads(); // wait for my block to finish
__threadfence(); // make sure global buffer writes are "visible"
// now do global sync
if (!threadIdx.x) { atomicAdd((int *)blkcnt, 1); } // mark my block done
if (!idx) { // am I the master block/thread?
while (*blkcnt < gridDim.x); // wait for all blocks to finish
*blkcnt = 0;
*bufrdy = BUFFER_FULL; // indicate that buffer is ready
__threadfence_system(); // push it out to mapped memory
completed_iter++;
}
}
}
int validate(const int *data, const int dsize, const int val) {
for (int i = 0; i < dsize; i++) {
if (data[i] != val) {
printf("mismatch at %d, was: %d, should be: %d\n", i, data[i], val);
return 0;
}
}
return 1;
}
int main() {
int *h_buf1, *d_buf1, *h_buf2, *d_buf2;
volatile int *m_bufrdy1, *m_bufrdy2;
// buffer and "signal" setup
cudaHostAlloc(&h_buf1, DSIZE * sizeof(int), cudaHostAllocDefault);
cudaHostAlloc(&h_buf2, DSIZE * sizeof(int), cudaHostAllocDefault);
cudaHostAlloc(&m_bufrdy1, sizeof(int), cudaHostAllocMapped);
cudaHostAlloc(&m_bufrdy2, sizeof(int), cudaHostAllocMapped);
cudaCheckErrors("cudaHostAlloc fail");
cudaMalloc(&d_buf1, DSIZE * sizeof(int));
cudaMalloc(&d_buf2, DSIZE * sizeof(int));
cudaCheckErrors("cudaMalloc fail");
cudaStream_t streamk, streamc;
cudaStreamCreate(&streamk);
cudaStreamCreate(&streamc);
cudaCheckErrors("cudaStreamCreate fail");
*m_bufrdy1 = BUFFER_EMPTY;
*m_bufrdy2 = BUFFER_EMPTY;
cudaMemset(d_buf1, 0xFF, DSIZE * sizeof(int));
cudaMemset(d_buf2, 0xFF, DSIZE * sizeof(int));
cudaCheckErrors("cudaMemset fail");
// inefficient crutch for choosing number of blocks
int nblock = 0;
cudaDeviceGetAttribute(&nblock, cudaDevAttrMultiProcessorCount, 0);
cudaCheckErrors("cudaDeviceGetAttribute fail");
test_kernel<<<nblock, nTPB, 0, streamk>>>(d_buf1, d_buf2, m_bufrdy1, m_bufrdy2, DSIZE);
cudaCheckErrors("kernel launch fail");
volatile int *bufrdy;
int *hbuf, *dbuf;
for (int i = 0; i < TOTAL_ITERS; i++) {
// ping pong on the host side
if (i & 1) {
bufrdy = m_bufrdy2;
hbuf = h_buf2;
dbuf = d_buf2;
} else {
bufrdy = m_bufrdy1;
hbuf = h_buf1;
dbuf = d_buf1;
}
while ((*bufrdy) != BUFFER_FULL);
cudaMemcpyAsync(hbuf, dbuf, DSIZE * sizeof(int), cudaMemcpyDeviceToHost, streamc);
cudaStreamSynchronize(streamc);
cudaCheckErrors("cudaMemcpyAsync fail");
*bufrdy = BUFFER_EMPTY; // release buffer back to device
if (!validate(hbuf, DSIZE, i)) {
printf("validation failure at iter %d\n", i);
exit(1);
}
}
printf("Completed %d iterations successfully\n", TOTAL_ITERS);
return 0;
}
上述代码中有几个值得注意的点：
指针传递
我们直接将cudaHostAlloc分配得到的指针传给了kernel，因为这块内存对于Host和GPU来说都是可见并且可直接在上面执行原子操作，并不需要额外的特殊处理。
stream并发
我们为kernel的执行和后面读回数据的memcpy创建了两个不同的stream，以保证他们之间不会产生相互阻塞。
block数量选择
在launch kernel的时候，我们直接启动了与设备上面SM个数相同的block，以保证这些block一定同时被调度并启动，不然全局同步可能会有问题。事实上，这种策略是最安全的但却不是最优的，因为多个block能够被调度到一个SM上面，这个上限并不是1。但是具体一个SM能够处理多少block非常取决于其资源占用，所以这个数字与kernel代码实现密切相关。此外，每个block内部固定启动256个线程，但是这个影响不大，因为反正这些线程是会被抢占调度的，但是block不会。
全局同步
在每个线程内部，处理完自己的数据以后，要先在block内部做一下所有线程同步，然后也要做一个全局所有block的同步来保证当前任务的所有数据都处理完成之后，才能去做下一项任务。由于CUDA编程模型的限制，这同样只能通过原子变量来实现。如果我们启动了数量大于设备最大能够调度个数的block，会导致这里在同步的时候并不是所有block都在运行，那么这里就会同步失败导致死循环。
同步开销
如果将workload设置为空，然后单独测量同步开销，我们会发现一些有意思的现象。整体上同步耗时与线程数量、block数量乘正相关。在只加载一个block，并且线程数<=64的时候，同步开销是最低的，这与我们对SM的理解一致：每个SM有两个wrap dispatcher，正好能够支持64个物理线程，如果再多就是串行执行了。另一方面，block数量与同步开销并不完全是正比关系。
生产者-消费者队列
（待完善）
一些思考
第一代Intel Xeon Phi计算加速卡（KNC）很神奇地采用了与GPU相似的工作模式，也就是用Host CPU来控制execution flow，让卡上那么多功能完备（性能是否跟CPU上一样就不好说了）的x86核心只负责参与计算，结果搞得那一代产品基本死掉。PS：天河2号上面全是这货（Xeon Phi 31S1P），目前来看基本上可以说是被坑了……
第二代KNL就聪明多了，直接就是一个众核的CPU，调度和计算都自己干，用户写程序的时候也不用为哪段代码跑在CPU哪段跑在加速卡上而操心，甚至MPI的已有程序都能无缝往上套，目前看起来还是取得了一定的市场认可。
所以个人猜想，是不是异构计算最终会趋于编程模型上的同构，即便处理器核是异构的，也能通过某些高性能的硬件互联机制，尽可能降低加速设备与逻辑控制单元之间的通信开销，从而让开发者用起来爽……
参考资料
- Nvidia Blog：Unified Virtual Addressing
- PT综述论文：Persistent-Threads-Style-Programming-Model-for-GPU-Computing，以及对应的简化版slides。
- 技术文档：Improving Real-Time Performance with CUDA Persistent Threads (CuPer) on the Jetson TX2
- 代码示例：A minimum CUDA persistent thread example
- 爆栈：host-device通过原子指令同步状态
- 爆栈：乒乓buffer实现kernel与host的数据交互