3.1 异构计算
在异构计算系统上进行的并行计算通常称为异构计算(Heterogeneous Computing)。它能够经济有效地实现高计算能力，可扩展性强，能够非常高效地利用计算资源。典型的例子就是Apple Inc.的异构芯片M1。
CUDA是一种实现异构计算的编程模型，开发过程中基本是围绕Host和Device展开的。
- 主机（Host）：CPU和内存（host memory），主要负责逻辑控制；
- 设备（Device）：GPU和显存（device memory），主要负责大数据计算。
在编程变量命令中常用‘h’和‘d’来区分内存（主存）和显存。
3.2 CUDA程序开发
3.2.1 CUDA安装
CUDA适用于所有包含NVIDIA GPU的服务器，工作站，个人电脑，嵌入式平台等电子设备。
安装方法：
- Linux：按照教程cuda-installation-guide-linux，需要6/7个步骤即可。
- Jetson：直接利用NVIDIA SDK MANAGER 或者SD image进行刷机即可，详见jetpack。
查看当前设备中的GPU状态：
- 服务器，工作站，个人电脑：
nvidia-smi； - Jetson等设备：
jtop
查看当前设备参数：执行deviceQuery程序。以Ubuntu为例，deviceQuery程序在：/usr/local/cuda/samples/1\_Utilities/deviceQuery。
3.2.2 CUDA程序编写
用C语言编写CUDA程序时，并行代码和串行代码是可以放在一起的。
编写思路：
- 把输入数据从CPU内存复制到GPU显存。
- 在执行芯片上缓存数据，加载GPU程序并执行。
- 将计算结果从GPU显存中复制到CPU内存中。
CUDA的C语言开发采用Extended C模式，首先介绍三个重要的关键字——执行空间说明符：
-
__global__执行空间说明符将函数声明为内核（Kernel），被其声明的函数：
- 在设备上执行；
- 可从主机调用，可在计算能力为3.2或更高的设备调用；
- 在调用时必须指定其执行配置；
- 是异步的，这意味着它在设备完成执行之前已经返回。
- 必须具有void返回类型，并且不能是类的成员。
-
__device__声明的函数：
- 在设备上执行；
- 只能从设备调用；
- 不能和
__global__一起使用。
-
__host__声明的函数：
- 在主机上执行；
- 只能从主机调用；
- 不能和
__global__一起使用，但是能够和__device__一起使用，此情况下该函数同时被主机和设备编译。
// 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：调用核函数，指定其执行配置（execution configuration，定义调用了多少个线程）
VecAdd<<<1, N>>>(A, B, C);
...
}3.2.3 nvcc编译套件
nvcc是CUDA程序的编译器，位于bin/目录中。NVCC是一个编译器的调用集合，它会去调用很多其他的编译工具，比如gcc、cicc、ptxas、fatbinary等等。
CUDA程序是包含主机代码和设备代码的统一源代码。NVCC编译器在编译过程中将两者区分开来，主机代码是用ANSI C编写的简单代码，由主机的标准C编译器进一步编译，编译后以一个普通的CPU进程的方式运行。设备代码用ANSI C扩展语言编写，用关键词来标识数据并行函数Kernel以及与之相关的数据结构来扩展代码。设备代码通常由nvcc进一步编译并在GPU上执行。
NVCC的编译过程分为离线编译和即时编译两部分组成，如下图所示。
nvcc在为GPU进行编译的过程中有虚拟架构的概念，可以通过设置参数，来根据不同架构的GPU进行编译。nvcc编译具有向后兼容（Backwards compatibility）}性，现实中的GPU架构版本应该不低于虚拟架构的版本，否则无法运行。例如：
nvcc x.cu \
--gpu-architecture=compute_50 \
--gpu-code=sm_50, sm_52
# 上面的一定不要比下面的低。3.2.4 利用NVPROF查看程序执行情况
通常用来分析CUDA程序性能的工具有很多，包括nvprof和Visual Profiler。
- Visual Profiler是一种图形分析工具，能显示应用程序中CPU和GPU活动的时间线，并自动分析识别潜在的优化机会。
- nvprof也是用来测试与优化CUDA程序性能的工具。而nvprof没有可视化的图形界面，但能从命令行收集、查看和分析数据。
nvprof是其中非常方便的分析工具，适合用来入门。然而，目前主流的CUDA驱动不再支持nvprof命令，不过我们仍可以在 NVIDIA Nsight Systems 中使用。NVIDIA Nsight 系统是一种性能分析工具，旨在帮助开发人员在CPU和GPU之间调整和扩展软件。
nsys nvprof -o out.nvvp a.exe
# 查看所有函数在执行过程中对资源的占用情况
nsys nvprof -o out.nvvp --print-api-trace a.exe
# 查看核函数在执行过程中对资源的占用情况
nsys nvprof -o out.nvvp --print-gpu-trace a.exe3.3 线程层次
- Thread：Sequential execution Unit
- 所有线程执行相同的核函数；
- 并行执行。
- Thread Block：a group of threads
- 执行在一个Streaming Multiprocessor(SM, 流多处理器)；
- 同一个Block中的线程可以协作。
- Thread Grid：a collection of thread blocks
- 在一个Grid当中的Block可以在多个SM中执行。
每一块Block中线程的个数以及每一个Grid中Block的个数均可以用类型为int或dim3的数据在核函数的执行配置中进行指定。上图展示了一个二维Block以及二维Grid的例子。另外，同一个SM可以同时处理位于不同Block当中的threads或warps，这样能够最大化资源的使用效率。
dim3 grid[3, 2, 1], block[4, 3, 1]Block当中的threads的个数一般需要最好为32的倍数，因为warp是GPU上最小的执行粒度（finest possible granularity of execution），否则会存在资源浪费。
对于不同的设备，Grid可以最多配置的Block的数量以及每个维度上的数量限制、每个Block最多可配置的Thread的数量以及每个维度上的数量限制可以通过deviceQuery程序进行查阅（下图：deviceQuery on Jetson AGX Orin Developer Kit）。在一个Block中的线程之间可通过SM上的共享内存（shared memory）进行沟通。一个Grid中的不同Block通常更有可能被分配到不同的SM上（实际编程中并不需要考虑这种分配方式，CUDA会自动分配）
Block好像也不是必须的，为什么同时具有线程和Block两种层次？为什么不是只有线程？这样增加了一个层级的抽象，也增加了复杂度是为了什么？
Block层级的存在是与GPU的硬件架构密不可分的。前面有提到，一个Grid中的Block往往被分配到不同的SM当中，往后GPU慢慢地发展，只需要增加物理的SM，性能就可以无限地扩展。
【硬件调度】：
- Grid：共享同样的kernel和Context
- Block：SM级别的调度单位；
- Threads/Warp：CUDA core级别的调度单位。
【资源和通信】：- Grid：共享同样的kernel和Context
- Block：同一个SM， 同一个SM memory
- Threads/Warp：允许一个warp中的thread读取其他thread的值。
如何合理地设置Block数量与Grid数量？
block_size：32的倍数grid_size：(N + block_size -1) / block_size (这样可以保证最后配置到的线程数量超过数据集的大小N)
内置变量：在CUDA内核开始执行时，会自动对一组特殊的只读寄存器进行赋值，而这些寄存器所存值会为每个thread提供其所在Block在Grid中的索引（blockIdx.[x y (z)]），以及其在Block中的索引（threadIdx.[x y (z)]）。对于二维或三维的grids和blocks，赋值过程是按行优先次序进行的。通过引入索引的不同，使得执行核函数的线程之间产生不同，通过这些差异来确定各自处理哪些数据。
那我们如何能够得到一个线程在所有的线程中的索引值？比如：我们申请了4个线程块，每个线程块有8个线程，那么我们就申请了32个线程，那么我需要找到第3个线程块（编号为2的block）里面的第6个线程（编号为5的thread）在所有线程中的索引值怎么办？
这时我们就需要blockDim和gridDim这两个内置变量：
gridDim表示一个grid中包含多少个blockblockDim表示一个block中包含多少个线程
也就是说，在上面的那个例子中，gridDim.x=4,blockDim.x=8那么，我们要找的第22个线程（编号为21）的唯一索引就应该是，index = blockIdx.x * blockDim.x + threadIdx.x.
上图并不是说一个CUDA core只能对应一个thread，Jetson Nano只有128个CUDA core并不代表它最多可以进行128个线程。实际上每个CUDA core在完成当前thread之后，可以被分配到完成其他thread。
3.4 实验：向量并行加法
接下来，我们通过完成一个向量加法的实例来实践一下，我们来实现的cpu代码如下:
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
void add(const double *x, const double *y, double *z, const int N)
{
for (int n = 0; n < N; ++n)
{
z[n] = x[n] + y[n];
}
}
void check(const double *z, const int N)
{
bool has_error = false;
for (int n = 0; n < N; ++n)
{
if (fabs(z[n] - 3) > (1.0e-10))
{
has_error = true;
}
}
printf("%s\n", has_error ? "Errors" : "Pass");
}
int main(void)
{
const int N = 100000000;
const int M = sizeof(double) * N;
double *x = (double*) malloc(M);
double *y = (double*) malloc(M);
double *z = (double*) malloc(M);
for (int n = 0; n < N; ++n)
{
x[n] = 1;
y[n] = 2;
}
add(x, y, z, N);
check(z, N);
free(x);
free(y);
free(z);
return 0;
}3.4.1 GPU Code Workflow
为了完成这个程序，我们先要将数据传输给GPU，并在GPU完成计算的时候，将数据从GPU中传输给CPU内存。
- Allocate GPU Memory；
cudaMalloc((void **)&d_x, size);
cudaMalloc((void **)&d_y, size);
cudaMalloc((void **)&d_z, size);- Copy data from CPU to GPU；
cudaMemcpy(d_x, h_x, size, cudaMemcpyHostToDevice);
cudaMemcpy(d_y, h_y, size, cudaMemcpyHostToDevice);- Invoke the CUDA Kernel；
const int block_size = 128;
const int grid_size = (N + block_size - 1) / block_size;
add<<<grid_size, block_size>>>(d_x, d_y, d_z, N);- Copy result from GPU to CPU；
cudaMemcpy(h_z, d_z, size, cudaMemcpyDeviceToHost);- Release GPU Memory；
cudaFree(d_x);
cudaFree(d_y);
cudaFree(d_z);3.4.2
向量并行加法过程：
核函数：
void __global__ add(const double *x, const double *y, double *z, int count)
{
const int n = blockDim.x * blockIdx.x + threadIdx.x;
if( n < count)
{
z[n] = x[n] + y[n];
}
}
