CUDA C++ Programming Guide

1. 引言

在同等价位和power envelope情况下，Graphics Processing Unit(GPU) 比CPU提供了更高的instruction throughput 和 memory bandwidth。其它计算设备如FPGA，也是energy efficient，但是不如GPU的program flexibility。

GPU和CPU的设计目标不同：

• 1）CPU设计为：擅于execute a sequence of operations, called a thread, as fast as possible。且可同时并行运行几十个线程。
• 2）GPU设计为：擅于并行运行几千个线程（amortizing the slower single-thread performance to achieve greater throughput）。

GPU擅长高度并行计算，因此，将更多的transistors用于data processing（如浮点计算），而不是data caching和flow control：

GPU可hide memory access latencies with computation, instead of relying on large data caches and complex flow control来avoid long memory access latencies，因为data caches和flow control都是expensive in terms of transistors。

通常，应用中包含了串行部分和并行部分，而系统通常也设计为GPU和CPU的混合体，以实现总体性能最优。具有高度并行性的应用可利用GPU的并行特征实现 比 CPU更优的性能。

2. CUDA®: A General-Purpose Parallel Computing Platform and Programming Model

2006年11月，NVIDIA®引入CUDA®: A General-Purpose Parallel Computing Platform and Programming Model——可使用C++等语言作为high-level programming language。
可采购相应的显卡，也可选择GPU云服务器。

multicore CPU和manycore GPU的出现意味着主流处理器芯片都是并行系统。软件开发的挑战在于：transparently scale其并行性，以充分利用新增的processor cores。CUDA设计为降低C语言开发者的学习曲线。

GPU core有3个关键抽象：

• 1）a hierarchy of thread groups
• 2）shared memories
• 3）barrier synchronization

这些抽象提供了细粒度的data parallelism和thread parallelism，嵌套在粗粒度的data parallelism 和 task parallelism 之内。这将指导开发者，将问题分解为粗粒度的子问题，然后独立并行的借助blocks of threads来来解决。每个子问题分解为细粒度，可cooperatively parallel by all threads within the block来解决。这种分解有助于保持language expressivity，支持threads to cooperate when solving each sub-problem，同时支持自动扩展。
即：

• 粗粒度的子问题：可独立并行运算，使用blocks of threads。
• 每个子问题内的细粒度：由单个block内的的threads parallel来解决。可schedule到GPU内的任意可用multiprocessors上，以任意顺序，并行或串行均可。

2.1 kernels

CUDA中的函数不同于C++函数，称为kernels。
当kernels被调用时，可executed N times in parallel by N different CUDA threads，这一点不同于传统的C++函数。

kernel以__global__来标记，CUDA threads的数量由<<<...>>> execution configuration来指定。每个运行该kernel的thread都有a unique thread ID，可kernel可通过内置变量来访问。

如下例，借助内置变量threadIdx，实现两个size为$N$的向量$A+B=C$。

// Kernel definition
__global__ void VecAdd(float* A, float* B, float* C)
{
 int i = threadIdx.x;
 C[i] = A[i] + B[i];
}
int main()
{
 ...
 // Kernel invocation with N threads
 VecAdd<<<1, N>>>(A, B, C);
 ...
}

2.2 thread hierarchy

threadIdx为3-component vector，因此可使用one-dimensional thread index，或 two-dimensional thread index，或three-dimensional thread index来标记threads，所形成的的1-dimensional、2-dimensional或3-dimensional block of threads 称为a thread block。从而可invoke computation across the elements in a domain such as a vector, matrix, or volume。

以size为N*N的矩阵相加为例A+B=C：

// Kernel definition
__global__ void MatAdd(float A[N][N], float B[N][N],
 float C[N][N])
{
 int i = threadIdx.x;
 int j = threadIdx.y;
 C[i][j] = A[i][j] + B[i][j];
}
int main()
{
 ...
 // Kernel invocation with one block of N * N * 1 threads
 int numBlocks = 1;
 dim3 threadsPerBlock(N, N);
 MatAdd<<<numBlocks, threadsPerBlock>>>(A, B, C);
 ...
}

每个block内的thread数量是有限的，因为一个block内的所有thread运行在同一processor core上，必须共享该core上有限的内存资源。当前GPU，threadsPerBlock最多为1024。
但是，一个kernel可运行在多个block之上，总的threads数量等于numBlocks * threadsPerBlock。因此numBlocks取决于待处理的数据量。
blocks可以1维、2维或3维grid of thread blocks来表示：

当numBlocks>1时，需借助内置变量blockIdx和blockDim：

// Kernel definition
__global__ void MatAdd(float A[N][N], float B[N][N],
float C[N][N])
{
 int i = blockIdx.x * blockDim.x + threadIdx.x;
 int j = blockIdx.y * blockDim.y + threadIdx.y;
 if (i < N && j < N)
 C[i][j] = A[i][j] + B[i][j];
}
int main()
{
 ...
 // Kernel invocation
 dim3 threadsPerBlock(16, 16); //每个thread block size为16x16（256 threads）
 dim3 numBlocks(N / threadsPerBlock.x, N / threadsPerBlock.y);
 MatAdd<<<numBlocks, threadsPerBlock>>>(A, B, C);
 ...
}

thread block之间是相互独立的，可 以 任意顺序串行或并行执行。
block内的threads之间可通过共享内存来共享数据，可通过同步执行来coordinate memory accesses。
可通过调用__syncthreads()内置函数来指定kernel中的synchronization points。__syncthreads()可作为a barrier，该block内的所有threads必须wait before any is allowed to proceed。【__syncthreads()的作用为：实现block内线程间的同步。除__syncthreads()之外，Cooperative Groups API提供了更丰富的thread-synchronization primitives。】

为了高效协作，shared memory为离每个processor core很近的low-latency memory（更像是L1 cache），且__syncthreads()应为轻量级的。

2.3 memory hierarchy

CUDA threads运行时可从多个memory space中访问数据：

• 1）每个thread有其 private local memory。
• 2）每个thread block有共享内存，对该block内的所有threads可见，这些thread block共享内存具有与该block相同的生命周期。
• 3）所有的threads可访问相同的global memory。【此外，所有threads还额外有2个read-only内存空间：constant memory 和 texture memory。global、constant和texture 内存可根据不同的应用场景进行选择。这三个内存空间为persistent across kernel launches by the same application。texture memory可为特定格式的数据提供不同的addressing modes以及data filtering】

2.4 heterogeneous programming

CUDA编程模式中设定了：

• CUDA threads运行在物理隔离的device上（作为coprocessor），host上运行C++程序。即kernel运行在GPU上，而其他C++程序运行在CPU上。
• host和device维护各自独立的DRAM内存空间，分别为host memory和device memory。程序可通过CUDA runtime来管理对kernel可见的 global、constant和texture memory space：如device memory allocation、device memory deallocation以及data transfer between host and device memory。

Unified Memory可提供managed memory来bridge host memory space和device memory space。
Managed memory is accessible from all CPUs and GPUs in the system as a single, coherent memory image with a common address space.

3. gridDim VS blockIdx VS blockDim VS threadIdx

• gridDim：为grid dimension
• blockIdx：为grid内包含的block index
• blockDim：为block dimension
• threadIdx：为block内包含的thread index

默认，a kernel对应1个grid，可将该grid看成是三维立方体，gird内每个元素为一个block。若将grid描述为dim3 grid(10,10,2);，则该grid内一共有10*10*2=200个block。
每个block对应为一个三维的thread立方体。

通常，仅使用block立方体和thread立方体的x维度，对应为一维数组，此时，blockDim.x对应为每个block的size，gridDim.x为block的总数。
若launch kernel的参数为：

dim3 block_dim(128,1,1);
dim3 grid_dim(10,1,1);
kernel<<<grid_dim,block_dim>>>(...);

对应的：

• threadIdx.x的取值范围为[0,128)
• blockIdx.x的取值范围为[0,10)
• blockDim.x的值为128
• gridDim.x的值为10

若计算tid = threadIdx.x + blockIdx.x *blockDim.x，则对应为 [0,128) + 128 * [0,9)，对应tid的取值范围为[0,1280)。当需要map threads to tasks时，这种方式可提供unique identifier for all of your threads in your kernel。

参考资料

[1] 《CUDA C++ Programming Guide》
[2] NVCC
[3] 《Matrix computations on the GPU CUBLAS, CUSOLVER and MAGMA by example》
[4] 《Flexible matrix multiplication kernels on GPUs》
[5] Cuda gridDim and blockDim