【参加CUDA线上训练营】CUDA并行计算基础 - Chapter 3- 初识CUDA

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），被其声明的函数：
• 1. 在设备上执行；
• 2. 可从主机调用，可在计算能力为3.2或更高的设备调用；
• 3. 在调用时必须指定其执行配置；
• 4. 是异步的，这意味着它在设备完成执行之前已经返回。
• 5. 必须具有void返回类型，并且不能是类的成员。
• __device__ 声明的函数：
• 1. 在设备上执行；
• 2. 只能从设备调用；
• 3. 不能和__global__一起使用。
• __host__ 声明的函数：
• 1. 在主机上执行；
• 2. 只能从主机调用；
• 3. 不能和__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.exe

3.3 线程层次

• Thread：Sequential execution Unit
• 1. 所有线程执行相同的核函数；
• 2. 并行执行。
• Thread Block：a group of threads
• 1. 执行在一个Streaming Multiprocessor(SM, 流多处理器)；
• 2. 同一个Block中的线程可以协作。
• Thread Grid：a collection of thread blocks
• 1. 在一个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中包含多少个block
• blockDim表示一个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];
    }

}