GPU 硬件与 CUDA 程序开发工具

GPU 硬件

在由 CPU 和 GPU 构成的异构计算平台中，通常将起控制作用的 CPU 称为 主机（host），
将起加速作用的 GPU 称为 设备（device）。

主机和设备都有自己的 DRAM，之间一般由 PCIe 总线连接。

GPU 计算能力不等价于计算性能；表征计算性能的一个重要参数是 浮点数运算峰值（FLOPS）。
浮点数运算峰值有单精度和双精度之分。对于 Tesla 系列的 GPU，双精度下 FLOPS 一般是单精度下的 1/2;
对于 GeForce 系列的 GPU，双精度下 FLOPS 一般是单精度下的 1/32。

影响计算性能的另一个参数是 GPU 内存带宽（显存）。

CUDA 程序开发工具

CUDA；
OpenCL，更为通用的各种异构平台编写并行程序的框架，AMD 的 GPU 程序开发工具；
OpenACC，由多公司共同开发的异构并行编程标准。

CUDA 提供两层 API，即 CUDA 驱动API 和 CUDA 运行时API。
CUDA 开发环境中，程序应用程序是以主机（CPU）为出发点的；应用程序可以调用 CUDA 运行时 API、
CUDA 驱动 API 和一些已有的 CUDA 库。

CUDA 开发环境搭建

linux 操作系统：linux下cuda环境搭建

windows10 操作系统：windows10下cuda环境搭建

nvidia-smi 检查与设置设备

nvidia-smi
    +-----------------------------------------------------------------------------+
    | NVIDIA-SMI 462.30       Driver Version: 462.30       CUDA Version: 11.2     |
    |-------------------------------+----------------------+----------------------+
    | GPU  Name            TCC/WDDM | Bus-Id        Disp.A | Volatile Uncorr. ECC |
    | Fan  Temp  Perf  Pwr:Usage/Cap|         Memory-Usage | GPU-Util  Compute M. |
    |                               |                      |               MIG M. |
    |===============================+======================+======================|
    |   0  GeForce MX450      WDDM  | 00000000:2B:00.0 Off |                  N/A |
    | N/A   39C    P8    N/A /  N/A |    119MiB /  2048MiB |      0%      Default |
    |                               |                      |                  N/A |
    +-------------------------------+----------------------+----------------------+
    
    +-----------------------------------------------------------------------------+
    | Processes:                                                                  |
    |  GPU   GI   CI        PID   Type   Process name                  GPU Memory |
    |        ID   ID                                                   Usage      |
    |=============================================================================|
    |  No running processes found                                                 |
    +-----------------------------------------------------------------------------+

CUDA Version， 11.2；
GPU Name，GeForce MX450，设备号为 0；如果系统中有多个 GPU 且只要使用其中某个特定的 GPU，
可以通过设置环境变量 CUDA_VISIBLE_DEVICES 的值，从而可以在运行 CUDA 程序前选定 GPU;
TCC/WDDM，WDDM（windows display driver model），其它包括 TCC（Tesla compute cluster）；
可以通过命令行 nvidia-smi -g GPU_ID -dm 0，设置为 WDDM 模式（1 为 TCC 模式）；
Compute mode, Default，此时同一个 GPU 中允许存在多个进程；其他模式包括 E.Process，
指的是独占进程模式，但不适用 WDDM 模式下的 GPU；
可以通过命令行 nvidia-smi -i GPU_ID -c 0，设置为 Default 模式（1 为 E.Process 模式）;
Perf，p8（GPU 性能状态，最大p0~最小p12）；

更多关于 nvidia-smi 的资料：nvidia-smi

CUDA 中的线程组织

CUDA 虽然支持 C++ 但支持得并不充分，导致 C++ 代码中有很多 C 代码的风格。

CUDA 采用 nvcc 作为编译器，支持 C++ 代码；nvcc 在编译 CUDA 程序时，
会将纯粹的 c++ 代码交给 c++ 编译器，自己负责编译剩下的 cu 代码。

C++ 的 Hello World 程序

1
2
3

>> g++ hello.cpp -o ./bin/hello.exe
>> ./bin/hello
msvc: hello world!

CUDA 的 Hello World 程序

使用 nvcc 编译纯粹 c++ 代码

1
2
3

>> nvcc -o ./bin/hello_cu.exe hello.cu 
>> ./bin/hello_cu.exe
nvcc: hello world!

在该程序中其实并未使用 GPU。

使用核函数的 CUDA 程序

一个利用了 GPU 的 CUDA 程序既有主机代码，又有设备代码（在设备中执行的代码）。
主机对设备的调用是通过 核函数（kernel function） 实现的。

int main()
{
    主机代码
    核函数的调用
    主机代码

    return 0；
}

核函数与 c++ 函数的区别：

必须加 __global__ 限定；
返回类型必须是空类型 void。

__global__ void hell_from__gpu()
{
    // 核函数不支持 c++ 的 iostream。
    printf("gpu: hello world!\n");
}

调用核函数的方式：

1	hello_from_gpu<<<1, 1>>>

主机在调用一个核函数时，必须指明在设备中指派多少线程。核函数中的线程常组织为若干线程块：

三括号中第一个数字是线程块的个数（number of thread block）；
三括号中第二个数字是每个线程块中的线程数（number of thread in per block）。

一个核函数的全部线程块构成一个网格（grid），线程块的个数称为网格大小（grid size）。每个线程块中含有相同数目的线程，该数目称为线程块大小（block size）。

所以，核函数的总的线程数即网格大小*线程块大小:

1	hello_from_gpu<<<grid size, block size>>>

调用核函数后，调用 CUDA 运行时 API 函数，同步主机和设备：

1	cudaDeviceSynchronize();

核函数中调用输出函数，输出流是先存放在缓冲区的，而缓冲区不会自动刷新。

CUDA 的线程组织

核函数的总线程数必须至少等于计算核心数时才有可能充分利用 GPU 的全部计算资源。

hello_from_gpu<<<2, 4>>>

网格大小是2，线程块大小是4，总线程数即8。核函数中代码的执行方式是 “单指令-多线程”，
即每个线程执行同一串代码。

从开普勒架构开始，最大允许的线程块大小是 2^10 (1024)，最大允许的网格大小是 2^31 - 1（一维网格）。

线程总数可以由两个参数确定：

gridDim.x, 即网格大小；
blockDim.x, 即线程块大小；

每个线程的身份可以由两个参数确定：

blockIdx.x, 即一个线程在一个网格中的线程块索引，[0, gridDm.x);
threadIdx.x, 即一个线程在一个线程块中的线程索引，[0, blockDim.x);

网格和线程块都可以拓展为三维结构（各轴默认为 1）：

三维网格 grid_size(gridDim.x, gridDim.y, gridDim.z);
三维线程块 block_size(blockDim.x, blockDim.y, blockDim.z);

相应的，每个线程的身份参数：

线程块ID (blockIdx.x, blockIdx.y, blockIdx.z);
线程ID (threadIdx.x, threadIdx.y, threadIdx.z);

多维网格线程在线程块上的 ID；

tid = threadIdx.z * (blockDim.x * blockDim.y)  // 当前线程块上前面的所有线程数
    + threadIdx.y * (blockDim.x)               // 当前线程块上当前面上前面行的所有线程数
    + threadIdx.x                              // 当前线程块上当前面上当前行的线程数

多维网格线程块在网格上的 ID:

1
2
3

bid = blockIdx.z * (gridDim.x * gridDim.y)
    + blockIdx.y * (gridDim.x)
    + blockIdx.x

一个线程块中的线程还可以细分为不同的 线程束（thread warp），即同一个线程块中
相邻的 warp_size 个线程（一般为 32）。

对于从开普勒架构到图灵架构的 GPU，网格大小在 x, y, z 方向的最大允许值为（2^31 - 1, 2^16 - 1, 2^16 -1）；
线程块大小在 x, y, z 方向的最大允许值为（1024， 1024， 64），同时要求一个线程块最多有 1024 个线程。

CUDA 的头文件

CUDA 头文件的后缀一般是 “.cuh”；同时，同时可以包含c/cpp 的头文件 “.h”、“.hpp”，采用 nvcc 编译器会自动包含必要的 cuda 头文件，如 , ，同时前者也包含了c++头文件。

使用 nvcc 编译 CUDA 程序

nvcc 会先将全部源代码分离为主机代码和设备代码；主机代码完整的支持 c++ 语法，而设备代码只部分支持。

nvcc 会先将设备代码编译为 PTX（parrallel thread execution）伪汇编代码，再将其编译为二进制 cubin目标代码。在编译为 PTX 代码时，需要选项 -arch=compute_XY 指定一个虚拟架构的计算能力；在编译为 cubin 代码时，需要选项 -code=sm_ZW 指定一个真实架构的计算能力，以确定可执行文件能够使用的 GPU。

真实架构的计算能力必须大于等于虚拟架构的计算能力，例如：

-arch=compute_35  -code=sm_60  (right)
-arch=compute_60  -code=sm_35  (wrong)

如果希望编译出来的文件能在更多的GPU上运行，则可以同时指定多组计算能力，例如：

-gencode arch=compute_35, code=sm_35
-gencode arch=compute_50, code=sm_50
-gencode arch=compute_60, code=sm_60

此时，编译出来的可执行文件将包含3个二进制版本，称为 胖二进制文件（fatbinary）。

同时，nvcc 有一种称为 实时编译（just-in-time compilation）机制，可以在运行可执行文件时从其中保留的PTX
代码中临时编译出一个 cubin 目标代码。因此，需要通过选项 -gencode arch=compute_XY, code=compute_XY，
指定所保留 PTX 代码的虚拟架构，例如：

-gencode arch=compute_35, code=sm_35
-gencode arch=compute_50, code=sm_50
-gencode arch=compute_60, code=sm_60  
-gencode arch=compute_70, code=compute_70

于此同时，nvcc 编译有一个简化的编译选项 -arch=sim_XY，其等价于：

-gencode arch=compute_XY, code=sm_XY  
-gencode arch=compute_XY, code=compute_XY

关于 nvcc 编译器的更多资料： nvcc。

简单 CUDA 程序的基本框架

单源文件 CUDA 程序基本框架

对于单源文件的 cuda 程序，基本框架为：

包含头文件

定义常量或宏

声明 c++ 自定义函数和 cuda 核函数的原型

int main()
{
    1. 分配主机和设备内存
    2. 初始化主机中数据
    3. 将某些数据从主机复制到设备
    4. 调用核函数在设备中计算
    5. 将某些数据从设备复制到主机
    6. 释放主机和设备内存
}

c++ 自定义函数和 cuda 核函数的定义

CUDA 核函数的要求：

返回类型必须是 void，但是函数中可以使用 return（但不可以返回任何值）；
必须使用限定符 __glolbal__，也可以加上 c++ 限定符；
核函数支持 c++ 的重载机制；
核函数不支持可变数量的参数列表，即参数个数必须确定；
一般情况下，传给核函数的数组（指针）必须指向设备内存（“统一内存编程机制”除外）；
核函数不可成为一个类的成员（一般以包装函数调用核函数，将包装函数定义为类成员）；
在计算能力3.5之前，核函数之间不能相互调用；之后，通过“动态并行”机制可以调用；
无论从主机调用还是从设备调用，核函数都在设备中执行（“<<<,>>>”指定执行配置）。

自定义设备函数

核函数可以调用不带执行配置的自定义函数，即 设备函数。

设备函数在设备中执行、在设备中被调用；而核函数在设备中执行、在主机中被调用。

__global__修饰的函数称为核函数，一般由主机调用、在设备中执行；
__device__修饰的函数称为设备函数，只能被核函数或其他设备函数调用、在设备中执行；
__host__修饰主机段的普通 c++ 函数，在主机中被调用、在主机中执行，一般可以省略；
可以同时用 __host__ 和 __device__ 修饰函数，从而减少代码冗余，此时编译器将
分别在主机和设备上编译该函数；
不能同时用 __global__ 和 __device__ 修饰函数；
不能同时用 __global__ 和 __host__ 修饰函数；
可以通过 __noinline__ 建议编译器不要将一个设备函数当作内联函数；
可以通过 __forceinline__ 建议编译器将一个设备函数当作内联函数。

设备函数可以有返回值。

CUDA 程序的错误检测

检测 CUDA 运行时错误的宏函数

定义检查 cuda 运行时 API 返回值 cudaError_t 的宏函数。

#define CHECK(call)                                                     \
do {                                                                    \
    const cudaError_t error_code = call;                                \
    if (error_code != cudaSuccess)                                      \
    {                                                                   \
        printf("CUDA ERROR: \n");                                       \
        printf("    FILE: %s\n", __FILE__);                             \
        printf("    LINE: %d\n", __LINE__);                             \
        printf("    ERROR CODE: %d\n", error_code);                     \
        printf("    ERROR TEXT: %s\n", cudaGetErrorString(error_code)); \
        exit(1);                                                        \
    }                                                                   \
}while(0);

因为核函数没有返回值，所以无法直接检查核函数错误。间接的方法是，在调用核函数后执行：

CHECK(cudaGetLastError());  // 捕捉同步前的最后一个错误。
CHECK(cudaDeviceSynchronize());  // 同步主机和设备。

核函数的调用是 异步的，即主机调用核函数后不会等待核函数执行完成、而是立刻执行之后的语句。同步操作较为耗时，一般尽量避免；同时，只要在核函数调用后还有对其他任何能返回错误值的 API 函数进行同步调用，都会触发主机和设备的同步并捕捉到核函数中可能发生的错误。

此外，主机和设备之间的数据拷贝会隐式地同步主机和设备。一般要获得精确的出错位置，还是需要显式地同步，例如调用 cudaDeviceSynchronize()。

或者，通过设置环境变量 CUDA_LAUNCH_BLOCKING 为 1，这样所有核函数的调用都将不再是异步的，而是同步的。就是说，主机调用一个核函数之后必须等待其执行完，才能向下执行。一般仅用于程序调试。

CUDA-MEMCHECK 检查内存错误

CUDA 提供了 CUDA-MEMCHECK 的工具集，包括 memcheck, racecheck, initcheck, synccheck.

>> cuda-memcheck --tool memcheck [options] app-name [options]

对于 memcheck 工具，可以简化为：

>> cuda-memcheck [options] app-name [options]

对于本例，可以通过如下方式检测错误：

    >> cuda-memcheck ./bin/check.exe
    ========= CUDA-MEMCHECK

CUDA ERROR: 
    FILE: check.cu
    LINE: 56
    ERROR CODE: 9
    ERROR TEXT: invalid configuration argument
    ========= Program hit cudaErrorInvalidConfiguration (error 9) due to "invalid configuration argument" on CUDA API call to     cudaLaunchKernel.
    =========     Saved host backtrace up to driver entry point at error
    =========     Host Frame:C:\Windows\system32\DriverStore\FileRepository\nvhq.inf_amd64_5550755be1247d27\nvcuda64.dll     (cuProfilerStop + 0x97c18) [0x2b8ca8]
    =========     Host Frame:C:\Windows\system32\DriverStore\FileRepository\nvhq.inf_amd64_5550755be1247d27\nvcuda64.dll     (cuProfilerStop + 0x9a2da) [0x2bb36a]
    =========     Host Frame:C:\Windows\system32\DriverStore\FileRepository\nvhq.inf_amd64_5550755be1247d27\nvcuda64.dll     [0x7b52e]
    =========     Host Frame:C:\Windows\system32\DriverStore\FileRepository\nvhq.inf_amd64_5550755be1247d27\nvcuda64.dll     (cuProfilerStop + 0x11ceaa) [0x33df3a]
    =========     Host Frame:C:\Windows\system32\DriverStore\FileRepository\nvhq.inf_amd64_5550755be1247d27\nvcuda64.dll     (cuProfilerStop + 0x137532) [0x3585c2]
    =========     Host Frame:D:\3_codes\CudaSteps\capter4\bin\check.exe [0x1679]
    =========     Host Frame:D:\3_codes\CudaSteps\capter4\bin\check.exe [0xd32b]
    =========     Host Frame:D:\3_codes\CudaSteps\capter4\bin\check.exe [0xd1a8]
    =========     Host Frame:D:\3_codes\CudaSteps\capter4\bin\check.exe [0xc6a1]
    =========     Host Frame:D:\3_codes\CudaSteps\capter4\bin\check.exe [0xcbf8]
    =========     Host Frame:D:\3_codes\CudaSteps\capter4\bin\check.exe [0xd944]
    =========     Host Frame:C:\Windows\System32\KERNEL32.DLL (BaseThreadInitThunk + 0x14) [0x17034]
    =========     Host Frame:C:\Windows\SYSTEM32\ntdll.dll (RtlUserThreadStart + 0x21) [0x52651]
    =========
    ========= Program hit cudaErrorInvalidConfiguration (error 9) due to "invalid configuration argument" on CUDA API call to     cudaGetLastError.
    =========     Saved host backtrace up to driver entry point at error
    =========     Host Frame:C:\Windows\system32\DriverStore\FileRepository\nvhq.inf_amd64_5550755be1247d27\nvcuda64.dll     (cuProfilerStop + 0x97c18) [0x2b8ca8]
    =========     Host Frame:C:\Windows\system32\DriverStore\FileRepository\nvhq.inf_amd64_5550755be1247d27\nvcuda64.dll     (cuProfilerStop + 0x9a2da) [0x2bb36a]
    =========     Host Frame:C:\Windows\system32\DriverStore\FileRepository\nvhq.inf_amd64_5550755be1247d27\nvcuda64.dll     [0x7b52e]
    =========     Host Frame:C:\Windows\system32\DriverStore\FileRepository\nvhq.inf_amd64_5550755be1247d27\nvcuda64.dll     (cuProfilerStop + 0x11ceaa) [0x33df3a]
    =========     Host Frame:C:\Windows\system32\DriverStore\FileRepository\nvhq.inf_amd64_5550755be1247d27\nvcuda64.dll     (cuProfilerStop + 0x137532) [0x3585c2]
    =========     Host Frame:D:\3_codes\CudaSteps\capter4\bin\check.exe [0x1461]
    =========     Host Frame:D:\3_codes\CudaSteps\capter4\bin\check.exe [0xcbfd]
    =========     Host Frame:D:\3_codes\CudaSteps\capter4\bin\check.exe [0xd944]
    =========     Host Frame:C:\Windows\System32\KERNEL32.DLL (BaseThreadInitThunk + 0x14) [0x17034]
    =========     Host Frame:C:\Windows\SYSTEM32\ntdll.dll (RtlUserThreadStart + 0x21) [0x52651]
    =========
    ========= ERROR SUMMARY: 2 errors

关于 CUDA-MEMCHECK 的更多内容，详见: CUDA-MEMCHECK。

获得 GPU 加速的关键

CUDA 事件计时

C++ 的计时方法：

GCC 和 MSVC 都有的 clock()函数；
原生的时间库；
GCC 的 gettimeofday()计时；
MSVC 的 QueryPerformanceCounter() 和 QueryPerformanceFrequency() 计时。

CUDA 基于 CUDA 事件的计时方法：

cudaEvent_t start, stop;
CHECK(cudaEventCreate(&start)); // 创建cuda 事件对象。
CHECK(cudaEventCreate(&stop));
CHECK(cudaEventRecord(start));  // 记录代表开始的事件。
cudaEventQuery(start);  // 强制刷新 cuda 执行流。

// run code.

CHECK(cudaEventRecord(stop));
CHECK(cudaEventSynchronize(stop)); // 强制同步，让主机等待cuda事件执行完毕。
float elapsed_time = 0;
CHECK(cudaEventElapsedTime(&curr_time, start, stop)); // 计算 start 和stop间的时间差（ms）。
printf("host memory malloc and copy: %f ms.\n", curr_time - elapsed_time);

由于 cuda 程序需要在主机和设备间传递数据，所以当计算强度较小时数据传输的性能对程序总耗时影响更大。
因此 cuda 的两种浮点数类型对程序性能的影响就较为明显。考虑提供编译选项，指定版本：

#ifdef USE_DP
    typedef double real;  // 双精度
    const real EPSILON = 1.0e-15;
#else
    typedef float real;   // 单精度
    const real EPSILON = 1.0e-6f;
#endif

在编译时，除了指定 GPU 计算能力 -arch=sm_50，还可以指定 c++ 优化等级 -O3;同时，可以指定其他
编译选项，如 -DUSE_DP 启用双精度版本。

>> nvcc -O3 -arch=sm_50 -DUSE_DP -o ./bin/clock.exe add.cu clock.cu main.cpp
...
>> ./bin/clock
using double precision version
host memory malloc and copy: 2.054112 ms.
device memory malloc: 9.063583 ms.
kernel function : 0.803360 ms.
cuda; no error
copy from device to host: 7.489505 ms.  

>> nvcc -O3 -arch=sm_50 -o ./bin/clock.exe add.cu clock.cu main.cpp
...
>> ./bin/clock     
host memory malloc and copy: 0.950240 ms.
device memory malloc: 5.298208 ms.
kernel function : 0.620512 ms.
cuda; no errors
copy from device to host: 3.034208 ms.

可见双精度版本基本上比单精度版本耗时多一倍。

nvprof 查看程序性能

>> nvprof ./bin/clock

如果没有输出结果，需要将nvprof的目录包含到环境环境变量中（不支持7.5 以上计算能力的显卡）。
推荐采用一代性能分析工具： Nvidia Nsight Systems.

影响 GPU 加速的关键因素

要获得可观的 GPU 加速，就必须尽量缩减主机和设备间数据传输所花时间的占比。

有些计算即使在 GPU 中速度不高也要尽量放在 GPU 中实现，以避免过多数据经由 PCIe 传递。

提高算术强度可以显著地提高 GPU 相对于 CPU 的加速比。

算术强度，是指一个计算问题中算术操作的工作量与必要的内存操作的工作量之比。对设备内存的访问速度取决于 GPU 的显存带宽。

核函数的并行规模。

并行规模可以用 GPU 中的线程数目来衡量。一个 GPU 由多个流多处理器SM（streaming multiprocessor）构成，每个 SM 中有若干 CUDA 核心。每个 SM 是相对独立的，一个 SM 中最多驻留的线程数一般为 2048 或 1024（图灵架构）。

若要 GPU 满负荷工作，则核函数中定义的线程总数要不少于某值，一般与 GPU 能够驻留的线程总数相当。

CUDA 的数学函数库

CUDA 提供的数学函数库提供了多种 数学函数，同时 CUDA 提供了一些高效率、低准确度的 内建函数。

CUDA 数学函数库的更多资料，详见：CUDA math.

CUDA 的内存组织

CUDA 中不同类型的内存

CUDA 中的内存类型有：全局内存、常量内存、纹理内存、寄存器、局部内存、共享内存。
CUDA 的内存，即设备内存，主机无法直接访问。

全局内存

全局内存（global memory），即核函数中所有线程都可以访问的内存，可读可写，由主机端分配和释放；如 cudaMalloc() 的设备内存 d_x, d_y, d_z。

全局内存由于没有放到 GPU 芯片上，所以具有较高的延迟和较低的访问速度，但是容量大（显存）。全局内存主要为核函数提供数据，并在主机和设备、设备和设备之间传递数据。

全局内存的生命周期由主机端维护，期间不同的核函数可以多次访问全局内存。

除以上动态分配的全局内存变量外，还可以使用 静态全局内存变量，其所占内存数量在编译器确定；这样的静态全局内存变量必须在所有主机和设备函数外部定义，例如：

1 2	__device__ real epsilon; // 单个静态全局内存变量， `__device` 表示是设备中的变量。 __device__ real arr[10]; // 固定长度的静态全局内存数组变量。

对于静态全局内存变量，其访问权限：

核函数中可以直接访问静态全局内存变量，不必以参数形式传给核函数；
主机中不可以直接访问静态全局内存变量，可以通过 cudaMemcpyToSymbol() 和 cudaMemcpyFromSymbol() 调用。

常量内存

常量内存（constant memory），仅有 64 kb，可见范围和生命周期与全局内存一样；具有缓存，从而高速；
常量内存仅可读、不可写。

使用常量内存的方法：一是在核函数外定义常量内存变量；二是向核函数传递常量参数，默认存放在常量内存：

核函数中可以直接访问常量全局内存变量，不必以参数形式传给核函数，但不可更改（只读）；
主机中不可以直接访问常量全局内存变量，可以通过 cudaMemcpyToSymbol() 和 cudaMemcpyFromSymbol() 调用。

纹理内存

纹理内存(texture memory)，类似常量内存，也是一种具有缓存的全局内存，具有相同可见范围和生命周期。

可以将某些只读的全局内存数据用 __ldg() 函数通过只读数据缓存（read-only data cache）读取，既可以达到使用纹理内存的加速效果，又可使代码简洁：

int __ldg(const int* ptr);  // 函数原型。

全局内存的读取在默认情况下就利用了 __ldg() 函数，所以不需要显式地使用。

寄存器

在核函数中定义的、不加任何限定符的变量一般存放在寄存器（register）；核函数中不加任何限定符的数组可能放在
寄存器，也可能放在局部内存中。寄存器可读可写。

各种内建变量，如 gridDim、blockDim 等都保存在特殊的寄存器中。

寄存器变量仅被一个线程看见，寄存器的生命周期也和所属线程相同。

寄存器内存在芯片上，是所有内存中访问速度最高的。一个寄存器占 32b（4字节），一个双精度浮点数占 2个寄存器。

局部内存

局部内存(local memory)也是全局内存的一部分，每个线程最多可以使用 512 kb 的局部内存，但过多使用会降低性能。
局部内存的用法类似寄存器。

共享内存

共享内存（shared memory）与寄存器类似，都是位于芯片上，读写速度较快。

共享内存对整个线程块可见，一个线程块上的所有线程都可以访问共享内存上的数据；共享内存的生命周期也与所属线程块一致。

共享内存的主要作用是减少对全局内存的访问，或者改善对全局内存的访问模式。

L1 和 L2 缓存

SM 层次的 L1 缓存（一级缓存）和设备层次 L2 缓存（二级缓存）。它们主要用来缓存全局内存和设备内存的访问。

SM 及其占有率

一个 GPU 由多个 SM（流多处理器）构成，一个 SM 包含如下资源：

一定数量的寄存器；
一定数量的共享内存；
常量内存的缓存；
纹理内存的缓存；
L1 缓存；
两个或四个线程束调度器，用于在不同线程上下文间迅速切换，及为准备就绪的线程束发出执行指令；
执行核心。

一般来说，要尽量让 SM 的占有率不小于某值（如 25%），才有可能获得较高的性能。

一个 SM 中最多拥有的线程块个数 Nb=16（开普勒和图灵架构）或 Nb=32（麦克斯韦、帕斯卡和伏特架构）；
一个 SM 中最多拥有的线程格式为 Nt=1028（图灵架构）或 Nt=2048（开普勒到伏特架构）。

在线程块中，每 32 个连续线程为一个 线程束。
SM 中线程的执行是以线程束为单位的，所以最好将线程块大小取为线程束大小（32个线程）的整数倍（如 128）.

CUDA 运行时 API 函数查询设备

使用 CUDA 运行时 API 函数查询所用GPU 规格。

#include "common/error.cuh"
#include <stdlib.h>


int main(int argc, char *argv[])
{
    int device_id = 0;
    if (argc > 1) device_id = atoi(argv[1]);

    CHECK(cudaSetDevice(device_id));

    cudaDeviceProp prop;
    CHECK(cudaGetDeviceProperties(&prop, device_id));

    printf("Device id: %d\n", device_id);
    printf("Device name: %s\n", prop.name);
    printf("Compute capability: %d.%d\n", prop.major, prop.minor);
    printf("Amount of global memory: %g GB\n", prop.totalGlobalMem/(1024.0*1024*1024));
    printf("Amount of constant memory: %g KB\n", prop.totalConstMem/1024.0);
    printf("Maximum grid size: %d, %d, %d\n", prop.maxGridSize[0], prop.maxGridSize[1], prop.maxGridSize[2]);
    printf("Maximum block size: %d, %d, %d\n", prop.maxThreadsDim[0],prop.maxThreadsDim[1],prop.maxThreadsDim[2]);
    printf("Number of SMs: %d\n", prop.multiProcessorCount);
    printf("Maximum amount of shared memory per block: %g KB\n", prop.sharedMemPerBlock/1024.0);
    printf("Maximum amount of shared memory per SM: %g KB\n", prop.sharedMemPerMultiprocessor/1024.0);
    printf("Maximum number of registers per block: %d K\n", prop.regsPerBlock/1024);
    printf("Maximum number of registers per SM: %d K\n", prop.regsPerMultiprocessor/1024);
    printf("Maximum number of threads per block: %d\n", prop.maxThreadsPerBlock);
    printf("Maximum number of threads per SM: %d\n", prop.maxThreadsPerMultiProcessor);

    return 0;
}

输出：

Device id: 0
Device name: GeForce MX450
Compute capability: 7.5
Amount of global memory: 2 GB
Amount of constant memory: 64 KB
Maximum grid size: 2147483647, 65535, 65535
Maximum block size: 1024, 1024, 64
Number of SMs: 14
Maximum amount of shared memory per block: 48 KB
Maximum amount of shared memory per SM: 64 KB
Maximum number of registers per block: 64 K
Maximum number of registers per SM: 64 K
Maximum number of threads per block: 1024
Maximum number of threads per SM: 1024

全局内存的合理使用

在各种设备内存中，全局内存具有最低的访问速度，往往是一个 CUDA 程序的性能瓶颈。

全局内存的合并与非合并访问

对全局内存的访问将触发内存事务，即数据传输。在启用了 L1 缓存的情况下，对全局内存的读取将首先尝试经过 L1 缓存；如果未命中，则尝试经过 L2 缓存；如果再次未命中，则直接从 DRAM 读取。

一次 数据传输处理 的数据量在默认情况下是 32 字节。
一次数据传输中，从全局内存转移到 L2 缓存的一片内存的首地址一定是 32 的整数倍。也就是说，一次数据传输只能从全局内存读取地址为 0-31 字节、32-63 字节等片段的数据。

合并度，即线程束请求的字节数与由此导致的所有内存事务中所传输的字节数之比。
如果所有数据传输处理的数据都是线程束所需要的，则合并度为 100%，即 合并访问；否则，即为 非合并访问。

以仅使用 L2 缓存的情况为例，一次数据传输指的就是将 32 字节数据从全局内存（DRAM）通过 32 字节的 L2 缓存片段（cache sector）传输到 SM。考虑一个线程束访问单精度浮点数类型的全局内存变量的场景，一个单精度浮点数占有 4 个字节，故一次访问需要 32*4 个字节的数据。在理想情况下，即合并度为 100% 时，将仅触发 128/32=4 次调用 L2 缓存的数据传输。如果线程束请求的全局内存地址刚好为 0-127 字节或 128-255 字节，就能与 4 次数据传输所处理的数据完全吻合，这种情况下就是合并访问。

64 位系统中基本数据类型的内存长度（字节）：

int size:               4
short size:             2
float size:             4
double size:            8
char size:              1
bool size:              1
long size:              4
int pointer size:       8
float pointer size:     8
double pointer size:    8
char pointer size:      8

矩阵转置

在核函数中，如果读取操作是非合并访问，则可以采用 只读数据缓存技术，通过加载函数 __ldg() 读取全局内存，从而对数据的读取进行缓存、缓解非合并访问的影响。

从帕斯卡架构开始，编译器会自动判断并调用 __ldg() 函数提升性能；对于开普勒架构、麦克斯韦架构，默认情况下不会使用 __ldg() 函数，需要手动配置。

对于核函数中全局内存的写入，则没有类似函数可用。所以若不能满足读取和写入都是合并的，一般应该尽量做到合并写入。

核函数中可以直接使用在函数外部由 #define 或 const 定义的常量，包括整型和浮点型常量。但是在windows平台下（MSVC编译器）核函数无法使用外部定义的 const 定义的浮点型常量。


#include "common/error.cuh"
#include <iostream>

#ifdef USE_DP
    typedef double real;  
    const real EPSILON = 1.0e-15;
#else
    typedef float real;   
    const real EPSILON = 1.0e-6f;
#endif

// using namespace std;  // 不能使用std，会导致 `copy()` 不能使用（命名冲突）。


__constant__ int TILE_DIM = 32;  // 设备内存中线程块中矩阵维度（线程块大小，最大1024）。

__global__ void copy(const real *src, real *dst, const int N);
__global__ void transpose1(const real *src, real *dst, const int N);
__global__ void transpose2(const real *src, real *dst, const int N);


int main()
{
    const int N = 10000; 
    const int M = N * N * sizeof(real);

    int SIZE = 0;
    CHECK(cudaMemcpyFromSymbol(&SIZE, TILE_DIM, sizeof(int)));   

    const int grid_size_x = (N + SIZE - 1)/SIZE; // 获取网格大小。
    const int grid_size_y = grid_size_x;

    const dim3 block_size(SIZE, SIZE);
    const dim3 grid_size(grid_size_x, grid_size_y);

    real *h_matrix_org, *h_matrix_res;
    h_matrix_org = new real[N*N];
    h_matrix_res = new real[N*N];
    for (int i = 0; i < N; ++i)
    {
        for (int j = 0; j < N; ++j)
        {
            h_matrix_org[j] = i;
        }   
    }

    float elapsed_time = 0;
    float curr_time = 0;
    cudaEvent_t start, stop;
    CHECK(cudaEventCreate(&start));
    CHECK(cudaEventCreate(&stop));
    CHECK(cudaEventRecord(start));
    cudaEventQuery(start);

    real *d_matrix_org, *d_matrix_res;
    CHECK(cudaMalloc(&d_matrix_org, M));
    CHECK(cudaMalloc(&d_matrix_res, M));
    CHECK(cudaMemcpy(d_matrix_org, h_matrix_org, M, cudaMemcpyDefault));

    copy<<<grid_size, block_size>>>(d_matrix_org, d_matrix_res, N);     
    CHECK(cudaMemcpy(h_matrix_res, d_matrix_res, M, cudaMemcpyDefault));

    CHECK(cudaEventRecord(stop));
    CHECK(cudaEventSynchronize(stop));
    CHECK(cudaEventElapsedTime(&curr_time, start, stop));    
    printf("matrix copy time cost: %f ms.\n", curr_time - elapsed_time);
    elapsed_time = curr_time;

    transpose1<<<grid_size, block_size>>>(d_matrix_org, d_matrix_res, N);     
    CHECK(cudaMemcpy(h_matrix_res, d_matrix_res, M, cudaMemcpyDefault));

    CHECK(cudaEventRecord(stop));
    CHECK(cudaEventSynchronize(stop));
    CHECK(cudaEventElapsedTime(&curr_time, start, stop));    
    printf("matrix transpose1 time cost: %f ms.\n", curr_time - elapsed_time);
    elapsed_time = curr_time;

    transpose2<<<grid_size, block_size>>>(d_matrix_org, d_matrix_res, N);     
    CHECK(cudaMemcpy(h_matrix_res, d_matrix_res, M, cudaMemcpyDefault));

    CHECK(cudaEventRecord(stop));
    CHECK(cudaEventSynchronize(stop));
    CHECK(cudaEventElapsedTime(&curr_time, start, stop));    
    printf("matrix transpose2 time cost: %f ms.\n", curr_time - elapsed_time);
    elapsed_time = curr_time;

    delete[] h_matrix_res;
    delete[] h_matrix_org;
    CHECK(cudaFree(d_matrix_org));
    CHECK(cudaFree(d_matrix_res));

    return 0;
}


__global__ void copy(const real *src, real *dst, const int N)
{
    // TILE_DIM = blockDim.x = blockDim.y
    const int nx = blockIdx.x * TILE_DIM + threadIdx.x; // 矩阵列索引。
    const int ny = blockIdx.y * TILE_DIM + threadIdx.y; // 矩阵行索引。
    const int index = ny * N + nx;

    if (nx >= N || ny >= N)
    {
        return;
    }

    dst[index] = src[index];  // 全局内存中数组也是线性存放的。
}

__global__ void transpose1(const real *src, real *dst, const int N)
{
    const int nx = threadIdx.x + blockIdx.x * TILE_DIM;
    const int ny = threadIdx.y + blockIdx.y * TILE_DIM;

    if (nx < N && ny < N)
    {
        // 矩阵转置（合并读取、非合并写入）。
        dst[nx*N + ny] = src[ny*N + nx];
    }
}

__global__ void transpose2(const real *src, real *dst, const int N)
{
    const int nx = threadIdx.x + blockIdx.x * TILE_DIM;
    const int ny = threadIdx.y + blockIdx.y * TILE_DIM;

    if (nx < N && ny < N)
    {
        // 矩阵转置（非合并读取、合并写入）。
        dst[ny*N + nx] = __ldg(&src[nx*N + ny]);   // 显示调用 `__ldg()` 函数缓存全局内存。 
    }
}

共享内存的合理使用

共享内存是一种可以被程序员直接操作的缓存，主要作用有两个：

减少核函数中对全局内存的访问次数，实现高效的线程块内部的通信；
提高全局内存访问的合并度。

数组归约

对于多线程程序，默认情况下不同线程的执行顺序是不固定的（线程间独立）。

采用 折半规约法，通过线程块对数据分片归约，最后再一并求和。

核函数中循环的每一轮都会被拆解、分配到线程块内的所有线程上执行，而不是一个线程连续执行一次完整循环。核函数中代码是 “单指令多线程” ，代码真正的执行顺序与出现顺序可能不同。所以线程 0、1、… 127之间实际上并行的。

保证一个线程块中所有线程在执行该语句后面的语句之前，都完全执行了前面的语句：通过 __syncthreads() 实现一个线程块中所有线程按照代码出现的顺序执行指令，但是不同线程块之间依然是独立、异步的。

共享内存变量，可以在核函数中通过限定符 __shared__ 定义一个共享内存变量，这样就相当于在每一个线程块中有一个该变量的副本。虽然每个副本都是独立的，但核函数中对共享变量的操作都将同时作用在所有副本上。

核函数中可以直接使用函数外部由 #define 或 const 定义的常量，但在 MSVC 中限制了核函数使用 const 定义的常量。

利用共享内存进行线程块之间的合作（通信）之前，都要进行同步，以确保共享内存变量中数据对于所有线程块内的所有线程都是准备好的。

共享内存的生命周期仅在核函数内，所以必须在核函数结束前将共享内存中需要的结果保存到全局内存。通过共享内存可以避免修改全局内存变量，同时不再要求全局内存数组为线程块大小的整数倍。

线程块的共享内存根据申请方式分为：静态共享内存变量和动态共享内存变量。前者在核函数中定义共享内存大小（通过编译期常量），后者在主机调用核函数时指定大小（可以提高可维护性）。

#include "../common/error.cuh"
#include "../common/floats.hpp"
#include <chrono>

using namespace std::chrono;

__constant__ int BLOCK_DIM = 128;


real reduce_cpu(const real *x, const int N)
{
    real sum = 0.0;
    for (int i = 0; i < N ; ++i)
    {
        sum += x[i];
    }

    return sum;
}

__global__ void reduce(real *x, real *y)
{
    // 这里执行迭代折半归约计算时，实际上的线程执行过程：
    // 1. 线程 0-127，offset = N/2, 迭代第一次；
    // 2. 线程 0-127，offset = N/4, 迭代第二次；
    // ...
    // 即，核函数中循环的每一轮都会被拆解、分配到线程块内的所有线程上执行，而不是一个  
    // 线程连续执行一次完整循环。
    const int tid = threadIdx.x;
    real *curr_x = x + blockIdx.x * blockDim.x;  // 当前线程块中处理的内存首地址。

    for (int offset = blockDim.x >> 1; offset > 0; offset >>=1) // 迭代折半归约。
    {
        // 由于条件筛选，实际导致每轮有效的线程数量减半，即 “线程束的分化”。
        // 要求数组大小为线程块大小的整数倍。
        if (tid < offset) 
        {
            // 核函数中代码是 “单指令多线程” ，代码真正的执行顺序与出现顺序可能不同。
            // 所以 线程 0、1、... 127之间实际上并行的。
            curr_x[tid] += curr_x[tid + offset];
        }

        // 保证一个线程块中所有线程在执行该语句后面的语句之前，都完全执行了前面的语句。
        // 实现一个线程块中所有线程按照代码出现的顺序执行指令，即线程1等待线程0，如此。
        // 但是不同线程块之间依然是独立、异步的。
        __syncthreads();
    }

    if (tid == 0)
    {
        // 通过线程块内同步，线程块 0 中的归约顺序：
        // 第一轮，curr_x[0] += curr_x[0+64], ... curr_x[63] += curr_x[63+64]；
        // 第二轮，curr_x[0] += curr_x[0+32], ... curr_x[31] += curr_x[31+32]；
        // 第三轮，curr_x[0] += curr_x[0+16], ... curr_x[15] += curr_x[15+16]； 
        // 第四轮，curr_x[0] += curr_x[0+ 8], ... curr_x[7 ] += curr_x[7 + 8]；  
        // 第五轮，curr_x[0] += curr_x[0+ 4], ... curr_x[3 ] += curr_x[3 + 4]； 
        // 第六轮，curr_x[0] += curr_x[0+ 2], curr_x[1 ] += curr_x[1 + 2]； 
        // 第七轮，curr_x[0] += curr_x[0+ 1]；           
        y[blockIdx.x] = curr_x[0];  
    }
}

__global__ void reduce_shared(const real *x, real *y, const int N)
{
    const int tid = threadIdx.x;
    const int bid = blockIdx.x;
    const int ind = bid * blockDim.x + tid;

    __shared__ real s_x[128];  // 定义线程块静态共享内存变量。 
    s_x[tid] = (ind < N) ? x[ind] : 0.0;  // 拷贝全局内存变量到线程块内的共享内存数据副本。
    __syncthreads();   // 同步线程块的数据拷贝操作，保证各线程块中数据对于块内线程都准备好。

    for (int offset = blockDim.x>>1; offset > 0; offset>>=1)
    {
        if (ind < offset)
        {
            s_x[tid] += s_x[tid + offset];
        }
        __syncthreads(); // 线程块内线程同步。
    }

    if (tid == 0)
    {
        y[bid] = s_x[0]; // 保存各个线程块中共享内存的0元素到全局内存。
    }
}

__global__ void reduce_shared2(const real *x, real *y, const int N)
{
    const int tid = threadIdx.x;
    const int bid = blockIdx.x;
    const int ind = bid * blockDim.x + tid;

    extern __shared__ real s_x[];  // 定义线程块动态共享内存变量，内存大小由主机调用核函数时定义。 
    s_x[tid] = (ind < N) ? x[ind] : 0.0;  // 拷贝全局内存变量到线程块内的共享内存数据副本。
    __syncthreads();   // 同步线程块的数据拷贝操作，保证各线程块中数据对于块内线程都准备好。

    for (int offset = blockDim.x>>1; offset > 0; offset>>=1)
    {
        if (ind < offset)
        {
            s_x[tid] += s_x[tid + offset];
        }
        __syncthreads(); // 线程块内线程同步。
    }

    if (tid == 0)
    {
        y[bid] = s_x[0]; // 保存各个线程块中共享内存的0元素到全局内存。
    }
}


int main()
{
    int N = 1e8;  // 单精度将发生 “大数吃小数” 的现象，导致结果完全错误；双精度没有问题。
    int M = N * sizeof(real);

    int block_size = 0;
    CHECK(cudaMemcpyFromSymbol(&block_size, BLOCK_DIM, sizeof(real)));
    int grid_size = (N + block_size - 1)/block_size; 

    real *h_x = new real[N];
    real *h_y = new real[grid_size]; 
    for (int i = 0; i < N; ++i)
    {
        h_x[i] = 1.23;
    }

    cout << FLOAT_PREC << endl;

    auto t1 = system_clock::now();

    // cpu归约，单精度下计算错误，大数吃小数。
    cout << "cpu reduce:  " << reduce_cpu(h_x, N) << endl;

    auto t2 = system_clock::now();
    double time = duration<double, std::milli>(t2 - t1).count();
    cout << "cpu reduce time cost: " << time << " ms" << endl;

    real *d_x, *d_y;
    int size = grid_size*sizeof(real);
    CHECK(cudaMalloc(&d_x, M)); 
    CHECK(cudaMalloc(&d_y, size)); // 数据分片后个线程块的归约结果数组。
    CHECK(cudaMemcpy(d_x, h_x, M, cudaMemcpyDefault)); 
    
    cudaEvent_t start, stop;
    CHECK(cudaEventCreate(&start));
    CHECK(cudaEventCreate(&stop));
    CHECK(cudaEventRecord(start));
    cudaEventQuery(start);

    // gpu归约，单精度下也能控制误差，稳健性更强。
    reduce<<<grid_size, block_size>>>(d_x, d_y); 
    CHECK(cudaMemcpy(h_y, d_y, size, cudaMemcpyDefault));
    CHECK(cudaGetLastError());

    float elap_time=0, curr_time=0;
    CHECK(cudaEventRecord(stop));
    CHECK(cudaEventSynchronize(stop));
    CHECK(cudaEventElapsedTime(&curr_time, start, stop));
    cout << "gpu reduce: " << reduce_cpu(h_y, grid_size) << endl;
    printf("gpu reduce time cost: %f ms\n", curr_time - elap_time);
    elap_time = curr_time;

    // gpu归约，采用静态共享内存的加速。
    reduce_shared<<<grid_size, block_size>>>(d_x, d_y, N); 
    CHECK(cudaMemcpy(h_y, d_y, size, cudaMemcpyDefault));
    CHECK(cudaGetLastError());

    CHECK(cudaEventRecord(stop));
    CHECK(cudaEventSynchronize(stop));
    CHECK(cudaEventElapsedTime(&curr_time, start, stop));
    cout << "gpu shared reduce: " << reduce_cpu(h_y, grid_size) << endl;
    printf("gpu shared reduce time cost: %f ms\n", curr_time - elap_time);
    elap_time = curr_time;    

    // gpu归约，采用动态共享内存的加速。
    // <<<grid_size, block_size, sharedMemSize>>>，第三个参数指定动态共享内存大小。
    int sharedMemSize = block_size * sizeof(real);  // 核函数中每个线程块的动态共享内存大小。
    reduce_shared2<<<grid_size, block_size, sharedMemSize>>>(d_x, d_y, N); 
    CHECK(cudaMemcpy(h_y, d_y, size, cudaMemcpyDefault));
    CHECK(cudaGetLastError());

    CHECK(cudaEventRecord(stop));
    CHECK(cudaEventSynchronize(stop));
    CHECK(cudaEventElapsedTime(&curr_time, start, stop));
    cout << "gpu shared2 reduce: " << reduce_cpu(h_y, grid_size) << endl;
    printf("gpu shared2 reduce time cost: %f ms\n", curr_time - elap_time);
    elap_time = curr_time; 

    delete[] h_x;
    delete[] h_y;
    CHECK(cudaFree(d_x));
    CHECK(cudaFree(d_y));

    return 0;
}

矩阵转置

由于共享内存访问速度快于全局内存，所以可以通过线程块内的共享内存将全局内存的非合并访问转为合并访问。

注意转置后的数组索引变换。

#include "../common/error.cuh"
#include "../common/floats.hpp"
#include <iomanip>
#include <string>
#include <fstream>

#define TILE_DIM  32

__constant__ int c_TILE_DIM = 32;  // 设备内存中线程块中矩阵维度（线程块大小，最大1024）。

void show(const real *matrix, const int N, std::string outfile, std::string title);
__global__ void transpose1(const real *src, real *dst, const int N);
__global__ void transpose2(const real *src, real *dst, const int N);
__global__ void transpose3(const real *src, real *dst, const int N);
__global__ void transpose4(const real *src, real *dst, const int N);



int main()
{
    // 由于显存 2 GB，float 为 4 字节，double 为 8 字节，所以在 transpose3, transpose4中：
    // float 矩阵维度不能超过 726；
    // double 矩阵维度不能超过 512；
    const int N = 500;  
    const int M = N * N * sizeof(real);

    int SIZE = 0;
    CHECK(cudaMemcpyFromSymbol(&SIZE, c_TILE_DIM, sizeof(int)));   

    const int grid_size_x = (N + SIZE - 1)/SIZE; // 获取网格大小。
    const int grid_size_y = grid_size_x;

    const dim3 block_size(SIZE, SIZE);
    const dim3 grid_size(grid_size_x, grid_size_y);

    real *h_matrix_org, *h_matrix_res;
    h_matrix_org = new real[N*N];
    h_matrix_res = new real[N*N];
    for (int i = 0; i < N; ++i)
    {
        for (int j = 0; j < N; ++j)
        {
            h_matrix_org[i * N + j] = i*1.0e-2;
        }   
    }
    // show(h_matrix_org, N, "result.txt", "origin matrix");
    
    real *d_matrix_org, *d_matrix_res;
    CHECK(cudaMalloc(&d_matrix_org, M));
    CHECK(cudaMalloc(&d_matrix_res, M));
    CHECK(cudaMemcpy(d_matrix_org, h_matrix_org, M, cudaMemcpyDefault));    

    float elapsed_time = 0;
    float curr_time = 0;
    cudaEvent_t start, stop;
    CHECK(cudaEventCreate(&start));
    CHECK(cudaEventCreate(&stop));
    CHECK(cudaEventRecord(start));
    cudaEventQuery(start);

    // 矩阵转置（全局内存合并读取、非合并写入）。
    transpose1<<<grid_size, block_size>>>(d_matrix_org, d_matrix_res, N);     
    CHECK(cudaMemcpy(h_matrix_res, d_matrix_res, M, cudaMemcpyDefault));
    // show(h_matrix_res, N, "result.txt", "transpose1");

    CHECK(cudaEventRecord(stop));
    CHECK(cudaEventSynchronize(stop));
    CHECK(cudaEventElapsedTime(&curr_time, start, stop));    
    printf("matrix transpose1 time cost: %f ms.\n", curr_time - elapsed_time);
    elapsed_time = curr_time;

    // 矩阵转置（全局内存非合并读取、合并写入）。
    transpose2<<<grid_size, block_size>>>(d_matrix_org, d_matrix_res, N);     
    CHECK(cudaMemcpy(h_matrix_res, d_matrix_res, M, cudaMemcpyDefault));
    // show(h_matrix_res, N, "matrix.txt", "transpose2");

    CHECK(cudaEventRecord(stop));
    CHECK(cudaEventSynchronize(stop));
    CHECK(cudaEventElapsedTime(&curr_time, start, stop));    
    printf("matrix transpose2 time cost: %f ms.\n", curr_time - elapsed_time);
    elapsed_time = curr_time;

    // 矩阵转置（通过共享内存全局内存合并读写）。
    transpose3<<<grid_size, block_size>>>(d_matrix_org, d_matrix_res, N);     
    CHECK(cudaMemcpy(h_matrix_res, d_matrix_res, M, cudaMemcpyDefault));
    // show(h_matrix_res, N, "result.txt", "transpose3");

    CHECK(cudaEventRecord(stop));
    CHECK(cudaEventSynchronize(stop));
    CHECK(cudaEventElapsedTime(&curr_time, start, stop));    
    printf("matrix transpose3 time cost: %f ms.\n", curr_time - elapsed_time);
    elapsed_time = curr_time;

    // 矩阵转置（通过共享内存、bank处理，实现全局内存合并读写）。
    transpose4<<<grid_size, block_size>>>(d_matrix_org, d_matrix_res, N);     
    CHECK(cudaMemcpy(h_matrix_res, d_matrix_res, M, cudaMemcpyDefault));
    // show(h_matrix_res, N, "result.txt", "transpose3");

    CHECK(cudaEventRecord(stop));
    CHECK(cudaEventSynchronize(stop));
    CHECK(cudaEventElapsedTime(&curr_time, start, stop));    
    printf("matrix transpose4 time cost: %f ms.\n", curr_time - elapsed_time);
    elapsed_time = curr_time;    

    delete[] h_matrix_res;
    delete[] h_matrix_org;
    CHECK(cudaFree(d_matrix_org));
    CHECK(cudaFree(d_matrix_res));

    return 0;
}


void show(const real *x, const int N, std::string outfile, std::string title)
{
    std::fstream out(outfile, std::ios::app);
    if (!out.is_open())
    {
        std::cerr << "invalid output file: " << outfile << endl;
        return;
    }

    out << "\n\n----------------" << title << endl;

    for (int i = 0; i < N; ++i)
    {
        out << endl;
        for (int j = 0; j < N; ++j)
        {
            out << std::setw(6) << x[i * N + j];
        }
    }
}

__global__ void transpose1(const real *src, real *dst, const int N)
{
    const int nx = threadIdx.x + blockIdx.x * c_TILE_DIM;
    const int ny = threadIdx.y + blockIdx.y * c_TILE_DIM;

    if (nx < N && ny < N)
    {
        // 矩阵转置（合并读取、非合并写入）。
        dst[nx*N + ny] = src[ny*N + nx];
    }
}

__global__ void transpose2(const real *src, real *dst, const int N)
{
    const int nx = threadIdx.x + blockIdx.x * c_TILE_DIM;
    const int ny = threadIdx.y + blockIdx.y * c_TILE_DIM;

    if (nx < N && ny < N)
    {
        // 矩阵转置（非合并读取、合并写入）。
        dst[ny*N + nx] = __ldg(&src[nx*N + ny]);   // 显示调用 `__ldg()` 函数缓存全局内存。 
    }
}

__global__ void transpose3(const real *src, real *dst, const int N)
{
    // 正常的做法中，全局内存的读写必有一个是非合并访问。
    // 现在通过将非合并访问转移到共享内存，利用共享内存的高性能（100倍全局内存），提高计算速度：  
    // 1. 首先将全局内存拷贝到线程块的共享内存；
    // 2. 然后从共享内存非合并访问，读取数据，合并写入全局内存。

    __shared__ real s_mat[TILE_DIM][TILE_DIM];  //二维静态共享内存，存储线程块内的一片矩阵。

    int bx = blockIdx.x * blockDim.x;  // 当前线程块首线程在网格中列索引。
    int by = blockIdx.y * blockDim.y;  // 当前线程块首线程在网格中行索引。

    int tx = threadIdx.x + bx;  // 当前线程在网格中列索引。
    int ty = threadIdx.y + by;  // 当前线程在网格中行索引。

    if (tx < N && ty < N)
    {
        // 全局内存合并访问，共享内存合并访问。
        s_mat[threadIdx.y][threadIdx.x] = src[ty * N + tx]; // 全局内存中二维矩阵一维存储。
    }
    __syncthreads();
    
    // 全局内存合并访问。
    if (tx < N && ty < N)
    {
        // 局部矩阵转置和全局内存合并写入。
        int x = by + threadIdx.x;
        int y = bx + threadIdx.y;
        dst[y * N + x] = s_mat[threadIdx.x][threadIdx.y]; 
    }
}

__global__ void transpose4(const real *src, real *dst, const int N)
{
    // 通过修改数组行大小，错开数组元素在共享内存bank中的分布，
    // 避免线程束的 32路bank冲突。
    __shared__ real s_mat[TILE_DIM][TILE_DIM + 1];  

    int bx = blockIdx.x * blockDim.x; 
    int by = blockIdx.y * blockDim.y; 

    int tx = threadIdx.x + bx;  
    int ty = threadIdx.y + by;  

    if (tx < N && ty < N)
    {
        s_mat[threadIdx.y][threadIdx.x] = src[ty * N + tx]; 
    }
    __syncthreads();
    
    if (tx < N && ty < N)
    {
        int x = by + threadIdx.x;
        int y = bx + threadIdx.y;
        dst[y * N + x] = s_mat[threadIdx.x][threadIdx.y]; 
    }
}

共享内存的 bank 冲突

共享内存在物理上被分为32个同样宽度（开普勒架构为 8 字节，其他为 4 字节）、能被同时访问的列向内存bank。

++++++++++++++++++++++++++++++++++++++

bank0     bank1    ...    bank31  

++++++++++++++++++++++++++++++++++++++

layer1    layer1   ...    layer1  
layer2    layer2   ...    layer2  
...  
layer32   layer32   ...   layer32

只要同一个线程束内的多个线程不同时访问同一个 bank 中不同层的数据，该线程束对共享内存的访问就只需要一次内存事务。当同一个线程束内的多个线程试图访问同一个 bank 中不同层的数据时，就会发生冲突。在同一线程束中的多个线程对同一个 bank 中的 n 层数据访问将导致 n 次内存事务，称为发生了 n 路 bank 冲突。

当线程束内的32个线程同时访问同一个 bank 的32个不同层，这将导致 32 路 bank 冲突。对于非开普勒架构，每个共享内存的宽带为 4 字节；于是每一层的32个 bank 将对应 32 个 float 数组元素。

使用共享内存来改善全局内存的访问方式不一定会提高核函数的性能；不要过早优化，在优化程序时要对不同的优化方案进行测试和比较。

原子函数的合理使用

cuda 中，一个线程的原子操作可以在不受其他线程的任何操作的影响下完成对某个（全局内存或共享内存）
数据的一套“读-改-写”操作。

完全在 GPU 中进行归约

有两种方法能够在GPU中得到最终结果：

用另一个核函数将较短的数组进一步归约；
在核函数末尾利用原子函数进行归约。

在代码实现中：

原子函数 atomicAdd(·)执行数组的一次完整的读-写操作；
传给 cudaMemcpy(·) 的主机内存可以是栈内存，也可以是堆内存；
主机函数可以和设备函数同名，但要遵循重载原则（参数列表不一致）。

#include "../common/error.cuh"
#include "../common/floats.hpp"
#include "../common/clock.cuh"


__global__ void reduce(real *x, real *y, const int N)
{
    int tid = threadIdx.x;
    int ind = tid + blockIdx.x * blockDim.x;

    extern __shared__ real curr_x[];
    curr_x[tid] = (ind < N) ? x[ind] : 0.0;    

    for (int offset = blockDim.x/2 ; offset > 0 ; offset /= 2)
    {
        if (tid < offset)
        {
            curr_x[tid] += curr_x[tid + offset];
        }
        __syncthreads();
    }

    if (tid == 0)
    {
        y[blockIdx.x] = curr_x[0];
    }
}

__global__ void reduce2(real *x, real *y, const int N)
{
    int tid = threadIdx.x;
    int ind = tid + blockIdx.x * blockDim.x;

    extern __shared__ real curr_x[];
    curr_x[tid] = (ind < N) ? x[ind] : 0.0;    

    for (int offset = blockDim.x/2 ; offset > 0 ; offset /= 2)
    {
        if (tid < offset)
        {
            curr_x[tid] += curr_x[tid + offset];
        }
        __syncthreads();
    }

    if (tid == 0)
    {
        // 原子函数 atomicAdd(*address, val).
        atomicAdd(y, curr_x[0]);
    }
}



int main()
{
    int N = 1e8;
    int M = N * sizeof(real);

    int bSize = 32;
    int gSize = (N + bSize - 1)/bSize;

    cout << FLOAT_PREC << endl;

    real *h_x, *h_y;
    h_x = new real[N];
    h_y = new real[gSize];
    for (int i = 0; i < N; ++i)
    {
        h_x[i] = 1.23;
    }

    cudaClockStart

    real *d_x, *d_y;
    CHECK(cudaMalloc(&d_x, M));
    CHECK(cudaMalloc(&d_y, gSize*sizeof(real)));
    CHECK(cudaMemcpy(d_x, h_x, M, cudaMemcpyDefault));

    cudaClockCurr

    reduce<<<gSize, bSize, (bSize+1)*sizeof(real)>>>(d_x, d_y, N);

    CHECK(cudaMemcpy(h_y, d_y, gSize*sizeof(real), cudaMemcpyDefault));
    real res = 0;
    for(int i = 0; i < gSize; ++i)
    {
        res += h_y[i];
    }
    cout << "reduce result: " << res << endl;

    cudaClockCurr

    reduce<<<gSize, bSize, (bSize)*sizeof(real)>>>(d_x, d_y, N);
    
    CHECK(cudaMemcpy(h_y, d_y, gSize*sizeof(real), cudaMemcpyDefault));
    res = 0.0;
    for(int i = 0; i < gSize; ++i)
    {
        res += h_y[i];
    }
    cout << "reduce result: " << res << endl;

    cudaClockCurr
    
    real *d_y2, *h_y2;
    h_y2 = new real(0.0);
    CHECK(cudaMalloc(&d_y2, sizeof(real)));

    // 采用原子函数、共享内存的核函数归约，
    // 由于减少了主机和设备间的数据传输，效率得以提高。
    reduce2<<<gSize, bSize, (bSize)*sizeof(real)>>>(d_x, d_y2, N);

    CHECK(cudaMemcpy(h_y2, d_y2, sizeof(real), cudaMemcpyDefault));
    cout << "reduce2 result: " << *h_y2 << endl;

    cudaClockCurr

    delete[] h_x;
    delete[] h_y;
    delete h_y2;
    CHECK(cudaFree(d_x));
    CHECK(cudaFree(d_y));
    CHECK(cudaFree(d_y2));

    return 0;
}

原子函数

原子函数对其第一个参数指向的数据进行一次“读-写-改”的原子操作，是不可分割的操作。第一个参数可以指向全局内存，也可以指向共享内存。

对所有参与的线程来说，原子操作是一个线程一个线程轮流进行的，没有明确的次序。原子函数没有同步功能。

原子函数的返回值为所指地址的旧值。

加法： T atomicAdd(T *address, T val)；
减法： T atomicSub(T *address, T val)；
交换： T atomicExch(T *address, T val);
最小值： T atomicMin(T *address, T val)；
最大值： T atomicMax(T *address, T val)；
自增： T atomicInc(T *address, T val)；
自减： T atomicDec(T *address, T val)；
比较交换： T atomicCAS(T *address, T compare, T val)；

邻居列表

两个粒子互为邻居的判断：他们的距离不大于一个给定的截断距离 rc。
基本算法：对每一个给定的粒子，通过比较它与所有其他粒子的距离来判断相应粒子对是否互为邻居。

#include "../common/error.cuh"
#include "../common/floats.hpp"
#include "../common/clock.cuh"
#include <fstream>
#include <regex>
#include <string>
#include <vector>


void read_data(const std::string &fstr, std::vector<real> &x, std::vector<real> &y);
void write_data(const std::string &fstr, const int *NL, const int N, const int M);
void find_neighbor(int *NN, int *NL, const real *x, const real *y, 
    const int N, const int M, 
    const real minDis);
__global__ void find_neighbor_gpu (int *NN, int *NL, const real *x, const real *y, 
    const int N, const int M, 
    const real mindDis);
__global__ void find_neighbor_atomic(int *NN, int *NL, const real *x, const real *y, 
    const int N, const int M, 
    const real minDis);


int main()
{
    cout << FLOAT_PREC << endl;

    std::string fstr = "xy.txt";
    std::string fout = "result.txt";
    std::vector<real> x, y;
    read_data(fstr, x, y);

    int N = x.size(), M = 10;
    real minDis = 1.9*1.9;

    int *NN = new int[N];
    int *NL = new int[N*M];
    for (int i = 0; i < N; ++i)
    {
        NN[i] = 0;
        for (int j = 0; j < M; ++j)
        {
            NL[i*M + j] = -1;
        }
    }

    int *d_NN, *d_NL;
    CHECK(cudaMalloc(&d_NN, N*sizeof(int)));
    CHECK(cudaMalloc(&d_NL, N*M*sizeof(int)));
    real *d_x, *d_y;
    CHECK(cudaMalloc(&d_x, N*sizeof(real)));
    CHECK(cudaMalloc(&d_y, N*sizeof(real)));

    cppClockStart

    find_neighbor(NN, NL, x.data(), y.data(), N, M, minDis);
    // write_data(fout, NL, N, M);
    cppClockCurr

    cudaClockStart

    CHECK(cudaMemcpy(d_x, x.data(), N*sizeof(real), cudaMemcpyDefault));
    CHECK(cudaMemcpy(d_y, y.data(), N*sizeof(real), cudaMemcpyDefault));

    int block_size = 128;
    int grid_size = (N + block_size - 1)/block_size;
    find_neighbor_atomic<<<grid_size, block_size>>>(d_NN, d_NL, d_x, d_y, N, M, minDis);

    CHECK(cudaMemcpy(NN, d_NN, N*sizeof(int), cudaMemcpyDefault));
    CHECK(cudaMemcpy(NL, d_NL, N*M*sizeof(int), cudaMemcpyDefault));
    // write_data(fout, NL, N, M);

    cudaClockCurr

    CHECK(cudaMemcpy(d_x, x.data(), N*sizeof(real), cudaMemcpyDefault));
    CHECK(cudaMemcpy(d_y, y.data(), N*sizeof(real), cudaMemcpyDefault));
    find_neighbor_gpu<<<grid_size, block_size>>>(d_NN, d_NL, d_x, d_y, N, M, minDis);
    CHECK(cudaMemcpy(NN, d_NN, N*sizeof(int), cudaMemcpyDefault));
    CHECK(cudaMemcpy(NL, d_NL, N*M*sizeof(int), cudaMemcpyDefault));

    cudaClockCurr    

    write_data(fout, NL, N, M);

    delete[] NN;
    delete[] NL;
    CHECK(cudaFree(d_NN));
    CHECK(cudaFree(d_NL));
    CHECK(cudaFree(d_x));
    CHECK(cudaFree(d_y));

    return 0;
}


void find_neighbor(int *NN, int *NL, const real *x, const real *y, 
    const int N, const int M, 
    const real minDis)
{
    for (int i = 0; i < N; ++i)
    {
        NN[i] = 0;
    }

    for (int i = 0; i < N; ++i)
    {
        for (int j = i + 1; j < N; ++j)
        {
            real dx = x[j] - x[i];
            real dy = y[j] - y[i];
            real dis = dx * dx + dy * dy;
            if (dis < minDis) // 比较平方，减少计算量。
            {
                NL[i*M + NN[i]] = j;  // 一维数组存放二维数据。
                NN[i] ++;
                NL[j*M + NN[j]] = i;  // 省去一般的判断。
                NN[j]++;
            }
        }
    }
}

__global__ void find_neighbor_gpu (int *NN, int *NL, const real *x, const real *y, 
    const int N, const int M, 
    const real minDis)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;

    if (i < N)
    {
        int count = 0;  // 寄存器变量，减少对全局变量NN的访问。
        for (int j = 0; j < N; ++j)  // 访问次数 N*N，性能降低。
        {
            real dx = x[j] - x[i];
            real dy = y[j] - y[i];
            real dis = dx * dx + dy * dy;

            if (dis < minDis && i != j) // 距离判断优先，提高“假”的命中率。
            {
                // 修改了全局内存NL的数据排列方式，实现合并访问（i 与 threadIdx.x的变化步调一致）。
                // ???
                NL[(count++) * N + i] = j;
            }
        }

        NN[i] = count;
    }
}

__global__ void find_neighbor_atomic(int *NN, int *NL, const real *x, const real *y, 
    const int N, const int M, 
    const real minDis)
{
    // 将 cpu 版本的第一层循环展开，一个线程对应一个原子操作。
    int i = blockIdx.x * blockDim.x + threadIdx.x;

    if (i < N)
    {
        NN[i] = 0;

        for (int j = i + 1; j < N; ++j)
        {
            real dx = x[j] - x[i];
            real dy = y[j] - y[i];
            real dis = dx * dx + dy*dy;
            if (dis < minDis)
            {
                // 原子函数提高的性能，但是在NL中产生了一定的随机性，不便于后期调试。
                int old_i_num = atomicAdd(&NN[i], 1);  // 返回值为旧值，当前线程对应点的邻居数
                NL[i*M + old_i_num] = j;  // 当前线程对应点的新邻居
                int old_j_num = atomicAdd(&NN[j], 1);  // 返回值为旧值，当前邻居点的邻居数
                NL[j*M + old_j_num] = i;  // 当前邻居点的新邻居
            }
        }
    }
}

void read_data(const std::string &fstr, std::vector<real> &x, std::vector<real> &y)
{
    x.clear();
    y.clear();

    std::fstream reader(fstr, std::ios::in);
    if (!reader.is_open())
    {
        std::cerr << "data file open failed.\n";
        return;
    }

    std::regex re{"[\\s,]+"};
    std::string line;
    while(std::getline(reader, line))
    {
       std::vector<std::string> arr{std::sregex_token_iterator(line.begin(), line.end(), re, -1), 
           std::sregex_token_iterator()};

		if (arr.size() < 2 || arr[0].find("#") != std::string::npos)
		{
			continue;
		}

		x.push_back(stod(arr[0]));
		y.push_back(stod(arr[1]));        
    }
}

void write_data(const std::string &fstr, const int *NL, const int N, const int M)
{
    std::fstream writer(fstr, std::ios::out);
    if (!writer.is_open())
    {
        std::cerr << "result file open failed.\n";
        return;
    }

    for (int i = 0; i < N; ++i)
    {
        writer << i << "\t";
        for (int j = 0; j < M; ++j)
        {
            int ind = NL[i*M + j];
            if (ind >= 0) 
            {
                writer << ind << "\t";
            }
        }

        writer << endl;
    }
}

线程束基本函数与协作组

线程束（warp），即一个线程块中连续32个线程。

单指令-多线程模式

一个GPU被分为若干个流多处理器（SM）。核函数中定义的线程块（block）在执行时将被分配到还没有完全占满的 SM。一个block不会被分配到不同的SM，同时一个 SM 中可以有多个 block。不同的block 之间可以并发也可以顺序执行，一般不能同步。当某些block完成计算任务后，对应的 SM 会部分或完全空闲，然后会有新的block被分配到空闲的SM。

一个 SM 以32个线程（warp）为单位产生、管理、调度、执行线程。
一个 SM 可以处理多个block，一个block可以分为若干个warp。

在同一时刻，一个warp中的线程只能执行一个共同的指令或者闲置，即单指令-多线程执行模型，（single instruction multiple thread, SIMT）。

当一个线程束中线程顺序的执行判断语句中的不同分支时，称为发生了 分支发散（branch divergence）。

if (condition)
{
    A;
}
else
{
    B;
}

首先，满足 condition 的线程或执行语句A，其他的线程会闲置；然后，不满足条件的将会执行语句B，其他线程闲置。当语句A和B的指令数差不多时，整个warp的执行效率就比没有分支的情况 低一半。

一般应当在核函数中尽量避免分支发散，但有时这也是不可避免的。如数组计算中常用的判断语句：

if(n < N)
{
    // do something.
}

该分支判断最多影响最后一个block中的某些warp发生分支发散，一般不会显著地影响性能。

有时能通过 合并判断语句 的方式减少分支发散；另外，如果两分支中有一个分支不包含指令，则即使发生分支发散也不会显著影响性能。

注意不同架构中的线程调度机制

线程束内的线程同步函数

__syncwarp(·)：当所涉及的线程都在一个线程束内时，可以将线程块同步函数 __syncthreads()
换成一个更加廉价的线程束同步函数__syncwarp(·)，简称 束内同步函数。

函数参数是一个代表掩码的无符号整型数，默认值是全部32个二进制位都为1，代表
线程束中的所有线程都参与同步。

关于 掩码(mask) 的简介文章：思否小姐姐：“奇怪的知识——位掩码”

协作组

协作组（cooperative groups），可以看作是线程块和线程束同步机制的推广，
提供包括线程块内部的同步与协作、线程块之间（网格级）的同步与协作、以及
设备与设备之间的同步与协作。

使用协作组需要包含如下头文件：

#include <cooperative_groups.h>
using namespace cooperative_groups;

线程块级别的协作组

协作组编程模型中最基本的类型是线程组 thread_group，其包含如下成员：

void sync()，同步组内所有线程；
unsigned size()，返回组内总的线程数目，即组的大小；
unsigned thread_rank()，返回当前调用该函数的线程在组内的标号（从0计数）；
bool is_valid()，如果定义的组违反了任何cuda限制，返回 false，否则true；

线程组类型有一个导出类型，线程块thread_block，其中定义了额外的函数：

dim3 group_index()，返回当前调用该函数的线程的线程块指标，等价于 blockIdx；
dim3 thread_index()，返回当前调用该函数的线程的线程指标，等价于 threadIdx；

通过 this_thread_block() 初始化一个线程块对象：

thread_block g = this_thread_block();  // 相当于一个线程块类型的常量。

此时，

g.sync() <===> __syncthreads()
g.group_index() <===> blockIdx
g.thread_index() <===> threadIdx

通过 tiled_partition() 可以将一个线程块划分为若干片（tile），每一片构成一个新的线程组。目前，仅支持将片的大小设置为 2 的整数次方且不大于 32。

thread_group g32 = tiled_partition(this_thread_block(), 32); // 将线程块划分为线程束。

可以继续将线程组划分为更细的线程组：

thread_group g4 = tiled_partition(g32, 4);

采用模板、在编译期划分 线程块片（thread block tile）：

thread_block_tile<32> g32 = tiled_partition<32>(this_thread_block());
thread_block_tile<32> g4 = tiled_partition<4>(this_thread_block());

线程块片具有额外的函数（类似线程束内函数）：

unsigned ballot(int predicate);
int all(int predicate);
int any(int predicate);
T shfl(T v, int srcLane);
T shfl_up(T v, unsigned d);
T shfl_down(T v, unsigned d);
T shfl_xor(T v, unsigned d);

与一般的线程束不同，线程组内的所有线程都要参与代码运行计算；同时，线程组内函数不需要指定宽度，因为该宽度就是线程块片的大小。

#include "../common/error.cuh"

const unsigned WIDTH = 8;
const unsigned BLOCK_SIZE = 16;
const unsigned FULL_MASK = 0xffffffff;

__global__ void test_warp_primitives(void)
{
    int tid = threadIdx.x;
    int laneId = tid % WIDTH;

    if (tid == 0)
    {
        printf("threadIdx.x: ");
    }
    printf("%2d  ", tid);
    if (tid == 0)
    {
        printf("\n");
    }

    if (tid == 0)
    {
        printf("laneId: ");
    }
    printf("%2d  ", laneId);
    if (tid == 0)
    {
        printf("\n");
    }

    // 从 FULL_MASK 出发， 计算 mask1（排除 0 号线程）掩码和mask2（仅保留 0 线程）掩码。
    unsigned mask1 = __ballot_sync(FULL_MASK, tid>0);
    unsigned mask2 = __ballot_sync(FULL_MASK, tid==0);

    if (tid == 0)
    {
        printf("FULL_MASK = %x\n", FULL_MASK);
    }
    if (tid == 1)
    {
        printf("mask1     = %x\n", mask1);
    }
    if (tid == 0)
    {
        printf("mask2     = %x\n", mask2);
    }

    // 从 FULL_MASK 出发计算线程束状态。
    // 因为不是所有线程的tid 都大于0，所以此处返回 0.
    int result = __all_sync(FULL_MASK, tid);
    if (tid == 0) 
    {
        printf("all_sync (FULL_MASK): %d\n", result);
    }
    // 从 mask1 出发计算线程束状态。
    // 因为mask1 中关闭了0号线程，所以剩下的所有线程tid > 0，此处返回 1.
    result = __all_sync(mask1, tid);
    if (tid == 1) 
    {
        printf("all_sync (mask1): %d\n", result);
    }

    // 从 FULL_MASK 出发计算线程束状态。
    // 因为存在线程的tid 都大于0，所以此处返回 1.
    result = __any_sync(FULL_MASK, tid);
    if (tid == 0)
    {
        printf("any_sync (FULL_MASK): %d\n", result);
    }
    // 从 mask2 出发计算线程束状态。
    // 因为mask2 中仅激活了 0 号线程，所以此处返回 0.
    result = __any_sync(mask2, tid);
    if (tid == 0)
    {
        printf("any_sync (mask2): %d\n", result);
    }

    // 从 FULL_MASK 出发，将每个线程束中 2号线程的tid广播到线程束内所有函数并作为返回值。
    // 所以在第一个线程束中，所有8个线程都将返回 laneId=2 线程（2 号线程）的tid值；
    // 在第二个线程束中，所有8个线程都也将返回 landId=2 线程（10 号线程）的tid值。
    int value = __shfl_sync(FULL_MASK, tid, 2, WIDTH);
    if (tid == 0)
    {
        printf("shfl:    ");
    }
    printf("%2d  ", value);
    if (tid == 0)
    {
        printf("\n");
    }

    // 从FULL_MASK出发，将每个线程束内 1-7 号线程取 0-6号线程的tid值并作为返回值。
    // 所以在第一个线程束中，0号线程返回自己的tid，1号线程返回0号线程的tid，2号线程返回1号线程tid, ...
    value = __shfl_up_sync(FULL_MASK, tid, 1, WIDTH);
    if (tid == 0)
    {
        printf("shfl_up:    ");
    }
    printf("%2d  ", value);
    if (tid == 0)
    {
        printf("\n");
    }

    // 从 FULL_MASK 出发，将每个线程束内 0-6号线程取 1-7 号线程的tid值并作为返回值。
    // 所以在第一个线程束中，0号线程返回1号线程的tid，2号线程返回3号线程的tid，..., 7号线程返回自己的tid。
    value = __shfl_down_sync(FULL_MASK, tid, 1, WIDTH);
    if (tid == 0)
    {
        printf("shfl_down:    ");
    }
    printf("%2d  ", value);
    if (tid == 0)
    {
        printf("\n");
    }    

    // 从 FULL_MASK 出发，将线程束中相邻的线程的tid相互传递并作为返回值。
    // 所以在第一个线程束中，0号线程返回1号线程的tid、1号线程返回0号线程的tid，2号线程返回3号线程的tid、  
    // 3号线程返回2号线程的tid，...
    value = __shfl_xor_sync(FULL_MASK, tid, 1, WIDTH);
    if (tid == 0)
    {
        printf("shfl_xor:    ");
    }
    printf("%2d  ", value);
    if (tid == 0)
    {
        printf("\n");
    }    
}


int main()
{
    test_warp_primitives<<<1, BLOCK_SIZE>>>();
    CHECK(cudaDeviceSynchronize());

    return 0;
}

数组归约程序的进一步优化

提高线程利用率

在当前的归约程序中，当 offset=64，只用了 1/2 的线程；当 offset=32，只用了 1/4 的线程；…
最终，当 offset=1，只用了 1/128 的线程；归约过程一共用了 log2(128) = 7 步，平均线程利用率 (1/2 + 1/4 + … + 1/128)/7 => 1/7。

而在归约前的数据拷贝中线程利用率为 100%，可以尽量把计算放在在归约前：让一个线程处理多个数据。

一个线程处理相邻若干个数据会导致全局内存的非合并访问。要保证全局内存的合并访问，这里需要保证相邻线程处理相邻数据，一个线程访问的数据需要有某种跨度。该跨度可以是线程块的大小，也可以是网格的大小；对于一维情况，分别是 blockDim.x 和 blockDim.x * gridDim.x。

避免反复分配与释放设备内存

设备内存的分配与释放是比较耗时的。通过采用静态全局内存替代动态全局内存，实现编译期的设备内存分配可以更加高效。

此外，应当尽量避免在较内存循环反复的分配和释放设备内存。

#include "../common/error.cuh"
#include "../common/floats.hpp"
#include "../common/clock.cuh"
#include <cooperative_groups.h>
using namespace cooperative_groups;

__constant__ unsigned FULL_MASK = 0xffffffff;
#define __gSize  10240
__device__ real static_y[__gSize];

__global__ void reduce_syncthreads(real *x, real *y, const int N);
__global__ void reduce_syncwarp(real *x, real *y, const int N);
__global__ void reduce_shfl_down(real *x, real *y, const int N);
__global__ void reduce_cp(real *x, real *y, const int N);
__global__ void reduce_cp_grid(const real *x, real *y, const int N);
real reduce_wrap(const real *x, const int N, const int gSize, const int bSize);
real reduce_wrap_static(const real *x, const int N, const int gSize, const int bSize);



int main()
{
    int N = 1e8;
    int M = N * sizeof(real);

    int bSize = 32;
    int gSize = (N + bSize - 1)/bSize;

    cout << FLOAT_PREC << endl;

    real *h_x, *h_x2, *h_y, *h_y2, *h_res;
    h_x = new real[N];
    h_x2 = new real[N];
    h_y = new real[gSize];
    h_y2 = new real[gSize];
    h_res = new real(0.0);
    for (int i = 0; i < N; ++i)
    {
        h_x[i] = 1.23;
        h_x2[i] = 1.23;
    }
    real initRes = 0.0;
    for (int i = 0; i < gSize ; ++i)
    {
        h_y2[i] = 0.0;
    }

    cudaClockStart

    real *d_x, *d_y, *d_res;
    CHECK(cudaMalloc(&d_x, M));
    CHECK(cudaMalloc(&d_y, gSize*sizeof(real)));
    CHECK(cudaMalloc(&d_res, sizeof(real)));
    CHECK(cudaMemcpy(d_x, h_x, M, cudaMemcpyDefault));

    cudaClockCurr

    reduce_syncthreads<<<gSize, bSize, (bSize)*sizeof(real)>>>(d_x, d_y, N);

    CHECK(cudaMemcpy(h_y, d_y, gSize*sizeof(real), cudaMemcpyDefault));
    real res = 0;
    for(int i = 0; i < gSize; ++i)
    {
        res += h_y[i];
    }
    cout << "reduce_syncthreads result: " << res << endl;
    cudaClockCurr
    
    CHECK(cudaMemcpy(d_res, &initRes, sizeof(real), cudaMemcpyDefault));
    reduce_syncwarp<<<gSize, bSize, bSize*sizeof(real)>>>(d_x, d_res, N);
    CHECK(cudaMemcpy(h_res, d_res, sizeof(real), cudaMemcpyDefault));
    cout << "reduce_syncwrap result: " << *h_res << endl;
    cudaClockCurr

    CHECK(cudaMemcpy(d_res, &initRes, sizeof(real), cudaMemcpyDefault));
    reduce_shfl_down<<<gSize, bSize, bSize*sizeof(real)>>>(d_x, d_res, N);
    CHECK(cudaMemcpy(h_res, d_res, sizeof(real), cudaMemcpyDefault));
    cout << "reduce_shfl_down result: " << *h_res << endl;
    cudaClockCurr

    CHECK(cudaMemcpy(d_res, &initRes, sizeof(real), cudaMemcpyDefault));
    reduce_cp<<<gSize, bSize, bSize*sizeof(real)>>>(d_x, d_res, N);
    CHECK(cudaMemcpy(h_res, d_res, sizeof(real), cudaMemcpyDefault));
    cout << "reduce_cp result: " << *h_res << endl;
    cudaClockCurr

    reduce_cp_grid<<<gSize, bSize, bSize*sizeof(real)>>>(d_x, d_y, N);
    CHECK(cudaMemcpy(h_y, d_y, gSize*sizeof(real), cudaMemcpyDefault));
    res = 0.0;
    for(int i = 0; i < gSize; ++i)
    {
        res += h_y[i];
    }    
    cout << "reduce_cp_grid result: " << res << endl;
    cudaClockCurr  

    res =  reduce_wrap(d_x, N, 10240, 128);
    cout << "reduce_wrap result: " << res << endl;
    cudaClockCurr  

    res =  reduce_wrap_static(d_x, N, 10240, 128);
    cout << "reduce_wrap_static result: " << res << endl;
    cudaClockCurr         

    delete[] h_x;
    delete[] h_y;
    delete h_res;
    CHECK(cudaFree(d_x));
    CHECK(cudaFree(d_y));
    CHECK(cudaFree(d_res));

    return 0;
}

__global__ void reduce_syncthreads(real *x, real *y, const int N)
{
    int tid = threadIdx.x;  // 线程块中线程在x方向的id。
    int ind = tid + blockIdx.x * blockDim.x; // 一维线程块中线程在GPU中的id。

    extern __shared__ real block_x[]; // 线程块共享内存。
    block_x[tid] = (ind < N)? x[ind] : 0;
    __syncthreads();  // 同步共享内存的拷贝操作，确保共享内存的数据已准备好。

    for(int offset = blockDim.x/2; offset > 0; offset /= 2)
    {
        if (tid < offset)
        {
            block_x[tid] += block_x[tid + offset];
        }
        __syncthreads(); // 同步线程块内线程。

    }    

    if (tid == 0)
    {
        y[blockIdx.x] = block_x[0];
    }

}

__global__ void reduce_syncwarp(real *x, real *y, const int N)
{
    const int tid = threadIdx.x;
    const int bid = blockIdx.x;
    const int ind = bid * blockDim.x + tid;

    extern __shared__ real block_arr[];
    block_arr[tid] = (ind < N) ? x[ind] : 0.0;
    __syncthreads();

    // 线程束之间的二分求和。
    for (int offset = blockDim.x/2; offset >= 32; offset /=2)  
    {
        if (tid < offset)
        {
            block_arr[tid] += block_arr[tid + offset];
        }
        __syncthreads(); // 同步线程块内的线程。
    }

    // 线程束内的二分求和。
    for (int offset = 16; offset > 0; offset /=2)
    {
        if (tid < offset)
        {
            block_arr[tid] += block_arr[tid + offset];
        }
        __syncwarp();  // 同步线程束内的线程。
    }

    if (tid == 0)
    {
        atomicAdd(y, block_arr[0]);  // 原子函数求和。
    }
}

__global__ void reduce_shfl_down(real *x, real *y, const int N)
{
    const int tid = threadIdx.x;
    const int bid = blockIdx.x;
    const int ind = bid * blockDim.x + tid;

    extern __shared__ real block_arr[];
    block_arr[tid] = (ind < N) ? x[ind] : 0.0;
    __syncthreads();

    for (int offset = blockDim.x /2 ; offset >= 32; offset /= 2)
    {
        if (tid < offset)
        {
            block_arr[tid] += block_arr[tid + offset];
        }

        __syncthreads();
    }

    // 在线程寄存器上定义一个变量y。
    real curr_y = block_arr[tid];

    for (int offset = 16; offset > 0; offset /= 2)
    {
        // 通过线程束洗牌函数，从FULL_MASK出发，
        // 将高线程号（数组索引）中的curr_y值平移到低线程号，通过设置偏移值为 offset，等价实现了线程束内的折半归约。
        curr_y += __shfl_down_sync(FULL_MASK, curr_y, offset);
    }  

    if (tid == 0)
    {
        atomicAdd(y, curr_y);
    }  
}

__global__ void reduce_cp(real *x, real *y, const int N)
{
    const int tid = threadIdx.x;
    const int bid = blockIdx.x;
    const int ind = bid * blockDim.x + tid;

    extern __shared__ real block_arr[];
    block_arr[tid] = (ind < N) ? x[ind] : 0.0;
    __syncthreads();

    for (int offset = blockDim.x /2 ; offset >= 32; offset /= 2)
    {
        if (tid < offset)
        {
            block_arr[tid] += block_arr[tid + offset];
        }

        __syncthreads();
    }

    real curr_y = block_arr[tid];   

    // 创建线程块片。
    thread_block_tile<32> g32 = tiled_partition<32>(this_thread_block());

    for (int offset = 16; offset > 0; offset /= 2)
    {
        // 线程块片的等价线程束内函数。
        curr_y += g32.shfl_down(curr_y, offset);
    }

    if (tid == 0)
    {
        atomicAdd(y, curr_y);
    }
}

__global__ void reduce_cp_grid(const real *x, real *y, const int N)
{
    const int tid = threadIdx.x;
    const int bid = blockIdx.x;
    extern __shared__ real block_arr[];

    real curr_y = 0.0;

    // 在归约前处理计算。
    // ???
    const int stride = blockDim.x * gridDim.x;
    for (int n = bid * blockDim.x + tid; n < N; n += stride)
    {
        curr_y += x[n];
    }

    block_arr[tid] = curr_y;
    __syncthreads();

    for (int offset = blockDim.x /2 ; offset >= 32; offset /= 2)
    {
        if (tid < offset)
        {
            block_arr[tid] += block_arr[tid + offset];
        }

        __syncthreads();
    }

    curr_y = block_arr[tid];
    thread_block_tile<32> g32 = tiled_partition<32>(this_thread_block());
    for (int offset = 16; offset > 0; offset /= 2)
    {
        // 线程块片的等价线程束内函数。
        curr_y += g32.shfl_down(curr_y, offset);
    }

    if (tid == 0)
    {
        y[bid] = curr_y;
    }    
}


real reduce_wrap(const real *x, const int N, const int gSize, const int bSize)
{
    const int ymem = gSize * sizeof(real);
    const int smem = bSize * sizeof(real);

    real h_y[1] = {0};
    real *d_y;
    CHECK(cudaMalloc(&d_y, ymem));

    // 使用两个核函数时，将数组 d_y 归约到最终结果的计算也是折半归约，  
    // 这比直接累加（使用原子函数或复制到主机再累加）要稳健（单精度下精度更高）。
    // 设备全局内存变量 d_x, d_y 对于每个线程块都是可见的，对于两个核函数是相同的。
    reduce_cp_grid<<<gSize, bSize, smem>>>(x, d_y, N);
    reduce_cp_grid<<<1, 1024, 1024*sizeof(real)>>>(d_y, d_y, gSize);

    CHECK(cudaMemcpy(h_y, d_y, sizeof(real), cudaMemcpyDefault));
    CHECK(cudaFree(d_y));

    return h_y[0];
}

real reduce_wrap_static(const real *x, const int N, const int gSize, const int bSize)
{
    real *d_y;
    CHECK(cudaGetSymbolAddress((void**)&d_y, static_y));   // 获取设备静态全局内存或常量内存的地址（指针）。

    reduce_cp_grid<<<gSize, bSize, bSize * sizeof(real)>>>(x, d_y, N);
    reduce_cp_grid<<<1, 1024, 1024*sizeof(real)>>>(d_y, d_y, gSize);

    real h_y[1] = {0};  
    CHECK(cudaMemcpy(h_y, d_y, sizeof(real), cudaMemcpyDefault));
    // CHECK(cudaFree(d_y));  // 全局内存由系统否则释放。

    return h_y[0];
}

CUDA 流

一个 CUDA 流一般是指由主机发出的、在设备中执行的cuda操作序列（即和cuda有关的操作，如主机—设备数据传输和核函数执行）。目前不考虑由设备段发出的流。

任何cuda操作都存在于某个cuda流，要么是 默认流（default stream），也称为空流；要么是明确指定的流。非默认的cuda流（非空流）都是在主机端产生与销毁。
一个cuda流由类型为 cudaStream_t 的变量表示，创建与销毁的方式：

cudaSteam_t stream;
CHECK(cudaStreamCreate(&stream));
...
CHECK(cudaStreamDestroy(stream));

主机中可以产生多个相互独立的cuda流，并实现cuda流之间的并行。

为了检查一个cuda流中所有操作是否都已在设备中执行完毕：

1 2	cudaError_t cudaStreamSynchronize(cudaStream_t stream); cudaError_t cudaStreamQuery(cudaStream_t stream);

cudaStreamSynchronize 会强制阻塞主机，直到其中的stream流执行完毕； cudaStreamQuery 不会阻塞主机，只是检查cuda流（stream）是否执行完毕，若是，则返回 cudaSuccess; 否则，返回 cudaErrorNotReady。

在默认流中重叠主机和设备计算

同一个cuda流在设备中都是顺序执行的。 在数组相加的例子中：

cudaMemcpy(d_x, h_x, M, cudaMemcpyDefault);
cudaMemcpy(d_y, h_y, M, cudaMemcpyDefault);
add<<<gridSize, blockSize>>>(d_x, d_y, d_z, N);
cudaMemcpy(h_z, d_z, M, cudaMemcpyDefault);

从设备的角度，以上4个cuda语句是按代码顺序执行的。

采用 cudaMemcpy 函数在主机与设备间拷贝数据，是具有隐式同步功能的。
所以从主机的角度看，数据传输是同步的或者说阻塞的，即主机在发出命令：

1	cudaMemcpy(d_x, h_x, M, cudaMemcpyDefault);

之后，会等待该命令执行完完毕，再接着往下走；数据传输时，主机是闲置的。
与此不同的是，核函数的启动是异步的或者说非阻塞的，即在主机发出命令：

1	add<<<gridSize, blockSize>>>(d_x, d_y, d_z, N);

之后，不会等待该命令执行完毕，而是立刻得到程序的控制权。紧接着发出：

1	cudaMemcpy(h_z, d_z, M, cudaMemcpyDefault);

然而，该命令不会被立刻执行，因为其与核函数同处默认流，需要顺序执行。

所以，主机在发出核函数调用后会立刻发出下一个命令；如果下一个命令是主机中的某个计算任务，那么主机就会在设备执行核函数的同时执行计算。这样就可以实现主机和设备的重叠计算。

当主机和设备的计算量相当时，将主机函数放在设备核函数后可以达到主机函数与设备函数并发执行的效果，从而有效地隐藏主机函数的执行时间。

非默认 cuda 流重叠多个核函数

要实现多个核函数之间的并行必须使用多个非默认 cuda 流。

使用多个流相对于使用一个流有加速效果；当流的数目超过某个阈值时，加速比就趋于饱和。制约加速比的因素：

GPU 计算资源，当核函数的线程总数超过某一值时，再增加流的数目就不会带来更高性能；
GPU 中能够并发执行的核函数的上限。

指定核函数的cuda流的方法：

1	kernal_func<<<grid_size, block_size, 0, stream>>>(params);

在调用核函数时，如果不需要使用共享内存，则该项设为0；同时指定cuda流的id。

计算能力为7，5的GPU能执行的核函数上限值为128。

非默认 cuda 流重叠核函数与数据传递

要实现核函数执行与数据传输的并发（重叠），必须让这两个操作处于不同的非默认流；同时，数据传输需要使用 cudaMemcpy 的异步版本 cudaMemcpyAsync。

异步传输由GPU的DMA（direct memory access）实现，不需要主机的参与。

使用异步的数据传输函数时，需要将主机内存定义为不可分页内存或者固定内存，从而防止在程序执行期间物理地址被修改。如果将可分页内存传递给 cudaMemcpyAsync 则会导致同步传输。

主机不可分页内存的分配与释放：

cudaError_t cudaMallocHost(void **ptr, size_t size);
或者
cudaError_t cudaHostAlloc(void **ptr, size_t size);

cudaError_t cudaFreeHost(void *ptr);

要利用多个流提升性能，一种方法是将数据和相应计算操作分为若干等分，然后在每个流中发布一个cuda操作序列。

如果核函数执行、主机与设备间的数据传输这3个cuda操作能完全并行执行，理论上最大加速比为 3。

#include "../common/error.cuh"
#include "../common/floats.hpp"
#include <math.h>
#include <stdio.h>


const int NUM_REPEATS = 10;
const int N1 = 1024;
const int MAX_NUM_STREAMS = 30;
const int N2 = N1 * MAX_NUM_STREAMS;
const int M2 = sizeof(real) * N2;
cudaStream_t streams[MAX_NUM_STREAMS];  // cuda流数组，全局变量由系统负责销毁。


const int N = 100000000;
const int M = sizeof(real) * N;
const int block_size = 128;
const int grid_size = (N - 1) / block_size + 1;


void timing(const real *h_x, const real *h_y, real *h_z,
    const real *d_x, const real *d_y, real *d_z,
    const int ratio, bool overlap);
void timing(const real *d_x, const real *d_y, real *d_z, 
    const int num);
void timing(const real *h_x, const real *h_y, real *h_z,
    real *d_x, real *d_y, real *d_z,
    const int num
);

int main(void)
{
    real *h_x = (real*) malloc(M);
    real *h_y = (real*) malloc(M);
    real *h_z = (real*) malloc(M);
    for (int n = 0; n < N; ++n)
    {
        h_x[n] = 1.23;
        h_y[n] = 2.34;
    }

    real *d_x, *d_y, *d_z;
    CHECK(cudaMalloc(&d_x, M));
    CHECK(cudaMalloc(&d_y, M));
    CHECK(cudaMalloc(&d_z, M));
    CHECK(cudaMemcpy(d_x, h_x, M, cudaMemcpyHostToDevice));
    CHECK(cudaMemcpy(d_y, h_y, M, cudaMemcpyHostToDevice));


    // host and kernal overlap.
    printf("Without CPU-GPU overlap (ratio = 10)\n");
    timing(h_x, h_y, h_z, d_x, d_y, d_z, 10, false);
    printf("With CPU-GPU overlap (ratio = 10)\n");
    timing(h_x, h_y, h_z, d_x, d_y, d_z, 10, true);

    printf("Without CPU-GPU overlap (ratio = 1)\n");
    timing(h_x, h_y, h_z, d_x, d_y, d_z, 1, false);
    printf("With CPU-GPU overlap (ratio = 1)\n");
    timing(h_x, h_y, h_z, d_x, d_y, d_z, 1, true);

    printf("Without CPU-GPU overlap (ratio = 1000)\n");
    timing(h_x, h_y, h_z, d_x, d_y, d_z, 1000, false);
    printf("With CPU-GPU overlap (ratio = 1000)\n");
    timing(h_x, h_y, h_z, d_x, d_y, d_z, 1000, true);


    // kernal and kernal overlap.
    for (int n = 0 ; n < MAX_NUM_STREAMS; ++n)
    {
        // 创建cuda流。
        CHECK(cudaStreamCreate(&(streams[n])));
    }

    for (int num = 1; num <= MAX_NUM_STREAMS; ++num)
    {
        timing(d_x, d_y, d_z, num);
    }

    for (int n = 0 ; n < MAX_NUM_STREAMS; ++n)
    {
        // 销毁cuda流。
        CHECK(cudaStreamDestroy(streams[n]));
    }


    // kernal and data transfering overlap.
    real *h_x2, *h_y2, *h_z2;
    CHECK(cudaMallocHost(&h_x2, M));
    CHECK(cudaMallocHost(&h_y2, M));
    CHECK(cudaMallocHost(&h_z2, M));
    for (int n = 0; n < N; ++n)
    {
        h_x2[n] = 1.23;
        h_y2[n] = 2.34;
    }

    for (int i = 0; i < MAX_NUM_STREAMS; i++)
    {
        CHECK(cudaStreamCreate(&(streams[i])));
    }

    for (int num = 1; num <= MAX_NUM_STREAMS; num *= 2)
    {
        timing(h_x2, h_y2, h_z2, d_x, d_y, d_z, num);
    }

    for (int i = 0 ; i < MAX_NUM_STREAMS; i++)
    {
        CHECK(cudaStreamDestroy(streams[i]));
    }

    CHECK(cudaFreeHost(h_x2));
    CHECK(cudaFreeHost(h_y2));
    CHECK(cudaFreeHost(h_z2));

    free(h_x);
    free(h_y);
    free(h_z);
    CHECK(cudaFree(d_x));
    CHECK(cudaFree(d_y));
    CHECK(cudaFree(d_z));

    return 0;
}

void cpu_sum(const real *x, const real *y, real *z, const int N_host)
{
    for (int n = 0; n < N_host; ++n)
    {
        z[n] = x[n] + y[n];
    }
}

void __global__ gpu_sum(const real *x, const real *y, real *z)
{
    const int n = blockDim.x * blockIdx.x + threadIdx.x;
    if (n < N)
    {
        z[n] = x[n] + y[n];
    }
}

void timing
(
    const real *h_x, const real *h_y, real *h_z,
    const real *d_x, const real *d_y, real *d_z,
    const int ratio, bool overlap
)
{
    float t_sum = 0;
    float t2_sum = 0;

    for (int repeat = 0; repeat <= NUM_REPEATS; ++repeat)
    {
        cudaEvent_t start, stop;
        CHECK(cudaEventCreate(&start));
        CHECK(cudaEventCreate(&stop));
        CHECK(cudaEventRecord(start));
        cudaEventQuery(start);

        if (!overlap)
        {
            cpu_sum(h_x, h_y, h_z, N / ratio);
        }

        gpu_sum<<<grid_size, block_size>>>(d_x, d_y, d_z);

        if (overlap)
        {
            // 主机函数与设备核函数重叠。
            cpu_sum(h_x, h_y, h_z, N / ratio);
        }
 
        CHECK(cudaEventRecord(stop));
        CHECK(cudaEventSynchronize(stop));
        float elapsed_time;
        CHECK(cudaEventElapsedTime(&elapsed_time, start, stop));
        printf("Time = %g ms.\n", elapsed_time);

        if (repeat > 0)
        {
            t_sum += elapsed_time;
            t2_sum += elapsed_time * elapsed_time;
        }

        CHECK(cudaEventDestroy(start));
        CHECK(cudaEventDestroy(stop));
    }

    const float t_ave = t_sum / NUM_REPEATS;
    const float t_err = sqrt(t2_sum / NUM_REPEATS - t_ave * t_ave);
    printf("Time = %g +- %g ms.\n", t_ave, t_err);
}

void __global__ add(const real *d_x, const real *d_y, real *d_z)
{
    const int n = blockDim.x * blockIdx.x + threadIdx.x;
    if (n < N1)
    {
        for (int i = 0; i < 100000; ++i)
        {
            d_z[n] = d_x[n] + d_y[n];
        }
    }
}

void timing(const real *d_x, const real *d_y, real *d_z, const int num)
{
    float t_sum = 0;
    float t2_sum = 0;

    for (int repeat = 0; repeat <= NUM_REPEATS; ++repeat)
    {
        cudaEvent_t start, stop;
        CHECK(cudaEventCreate(&start));
        CHECK(cudaEventCreate(&stop));
        CHECK(cudaEventRecord(start));
        cudaEventQuery(start);

        for (int n = 0; n < num; ++n)
        {
            int offset = n * N1;

            // 指定各个核函数的cuda流，实现核函数的并行。
            add<<<grid_size, block_size, 0, streams[n]>>>(d_x + offset, d_y + offset, d_z + offset);
        }
 
        CHECK(cudaEventRecord(stop));
        CHECK(cudaEventSynchronize(stop));
        float elapsed_time;
        CHECK(cudaEventElapsedTime(&elapsed_time, start, stop));

        if (repeat > 0)
        {
            t_sum += elapsed_time;
            t2_sum += elapsed_time * elapsed_time;
        }

        CHECK(cudaEventDestroy(start));
        CHECK(cudaEventDestroy(stop));
    }

    const float t_ave = t_sum / NUM_REPEATS;
    const float t_err = sqrt(t2_sum / NUM_REPEATS - t_ave * t_ave);
    printf("%g\n", t_ave);
}

void __global__ add2(const real *x, const real *y, real *z, int N)
{
    const int n = blockDim.x * blockIdx.x + threadIdx.x;
    if (n < N)
    {
        for (int i = 0; i < 40; ++i)
        {
            z[n] = x[n] + y[n];
        }
    }
}

void timing
(
    const real *h_x, const real *h_y, real *h_z,
    real *d_x, real *d_y, real *d_z,
    const int num
)
{
    int N1 = N / num;
    int M1 = M / num;

    float t_sum = 0;
    float t2_sum = 0;

    for (int repeat = 0; repeat <= NUM_REPEATS; ++repeat)
    {
        cudaEvent_t start, stop;
        CHECK(cudaEventCreate(&start));
        CHECK(cudaEventCreate(&stop));
        CHECK(cudaEventRecord(start));
        cudaEventQuery(start);

        for (int i = 0; i < num; i++)
        {
            int offset = i * N1;

            // 划分主机不可分页内存，实现异步的数据传输。
            // 每个cuda流都有各自的数据传输操作。
            CHECK(cudaMemcpyAsync(d_x + offset, h_x + offset, M1, 
                cudaMemcpyHostToDevice, streams[i]));
            CHECK(cudaMemcpyAsync(d_y + offset, h_y + offset, M1, 
                cudaMemcpyHostToDevice, streams[i]));

            int block_size = 128;
            int grid_size = (N1 - 1) / block_size + 1;

            // 指定核函数的cuda流。
            add2<<<grid_size, block_size, 0, streams[i]>>>(d_x + offset, d_y + offset, d_z + offset, N1);

            CHECK(cudaMemcpyAsync(h_z + offset, d_z + offset, M1, 
                cudaMemcpyDeviceToHost, streams[i]));
        }

        CHECK(cudaEventRecord(stop));
        CHECK(cudaEventSynchronize(stop));
        float elapsed_time;
        CHECK(cudaEventElapsedTime(&elapsed_time, start, stop));

        if (repeat > 0)
        {
            t_sum += elapsed_time;
            t2_sum += elapsed_time * elapsed_time;
        }

        CHECK(cudaEventDestroy(start));
        CHECK(cudaEventDestroy(stop));
    }

    const float t_ave = t_sum / NUM_REPEATS;
    const float t_err = sqrt(t2_sum / NUM_REPEATS - t_ave * t_ave);
    printf("%d %g\n", num, t_ave);
}

分子动力学模型

将静态函数放在头文件中，则该函数就有可能被编译为内联函数，从而提高效率。
适用于于需要被多个编译单元反复调用的函数。开发cuda程序时，也应该尽量优化对应的c++程序。
半步长推进。
粒子在t+dt时刻的坐标仅依赖t时刻的坐标、速度和力；但是t+dt时刻的速度依赖 t时刻的坐标、速度和t+dt时刻的力。所以首先，以t时刻的状态计算t+dt/2时刻的速度；然后计算t+dt时刻的坐标，同时更新t时刻的力到t+dt时刻；最后，以t+dt/2时刻的速度和t+dt时刻的力计算t+dt 时刻的速度。
常量内存比全局内存高速。
如果数据量在编译期就确定且不大（明显小于4KB），在核函数中仅被读取，而且一个线程束中的所有线程在某个时刻访问同一个地址，则该数据适合用传参的方式使用常量内存。
逐步分析程序的性能瓶颈，逐步优化。
将一个c++程序用cuda加速时，一般首先确定其中最耗时的部分并将其用cuda加速，从而快速提高程序性能。要得到最好的加速效果，需要尽可能多的将程序中可并行的计算用cuda加速。当然，在大数情况下，我们需要在付出和收获间找到一个平衡点。

cuda 标准库的使用

Thrust：类型 c++ 的标准模板库；
cuBLAS：基本线性代数子程序；
cuFFT：快速傅里叶变换；
cuSPARSE：稀疏矩阵；
cuRAND：随机数生成器；
cuSolver：稠密矩阵和稀疏矩阵计算库；
cuDNN：深度神经网络。

Thrust

Thrust：一个实现了众多基本并行算法的c++模板库，类似c++的标准库stl。

Thrust官方资料

数据结构

Thrust 中的数据结构主要是矢量容器，类似 stl 中的 std::vector:

存储于主机的容器 thrust::host_vector<typename>;
存储于设备的容器 thrust::device__vector<typename>;

容器的使用也类似于 stl：

// 包含头文件
#include<thrust/host_vector.h>
#include<thrust/device_vector.h>

// 定义并初始化主机内存
thrust::host_vector<double> arr(12, 0.0);

Thrust 函数可以直接调用设备上的矢量容器。

算法

Thrust 提供5类常用算法：变换，归约，前缀和，排序于搜索，选择性复制、替换、移除、分区等重排操作。

Thrust 函数的参数必须都来自于主机容器，或者都来自于设备容器；thrust::copy 除外。

如果程序中大量使用了 thrust 库，使用设备矢量较为合适；如果只是偶尔使用 Thrust 库，
则使用设备内存指针更为合适。

cuBLAS

cuBLAS，一个基本线性代数子程序，提供三层功能函数：

处理矢量之间的计算，如矢量之间的内积；
处理矩阵和矢量之间的运算，如相乘；
处理矩阵之间的运算，如相乘。

CUBLAS 中矩阵采用 列主序，即矩阵数据按列存储。

cuBLAS官方资料

cuSolver

cuSolver：稠密矩阵和稀疏矩阵计算库。

cuSolver 相比于 cuBLAS，专注于一些比较高级的线性代数计算，并由3个子库组成：

cuSolverDN，处理稠密矩阵线性代数计算；
cuSolverSP，处理稀疏矩阵线性代数计算；
cuSolverRF，处理稀疏矩阵分解。

cuSolver 库函数倾向于使用异步执行。为例保证一个 cuSolver 函数的工作已完成，
可以使用 cudaDeviceSynchronize() 函数进行同步。

cuSolver 中矩阵同样采用 列主序。

cuSolver官方资料

cuRAND

cuRAND：随机数生成器。

cuRAND 库中提供两种 API：主机API 和设备API。以主机API为例，使用方式：

1
2
3

#include <curand.h>

编译时指定链接选项 `-lcurand`

同时，主机API 分为两种使用方式：

使用设备产生伪随机数并存于设备数组；
使用主机产生伪随机数并存于主机数组。

cuRAND官方资料

GPU 硬件与 CUDA 程序开发工具

GPU 硬件

CUDA 程序开发工具

CUDA 开发环境搭建

nvidia-smi 检查与设置设备

CUDA 中的线程组织

C++ 的 Hello World 程序

123>> g++ hello.cpp -o ./bin/hello.exe>> ./bin/hellomsvc: hello world!

CUDA 的 Hello World 程序

使用 nvcc 编译纯粹 c++ 代码

使用 核函数 的 CUDA 程序

CUDA 的线程组织

CUDA 的头文件

使用 nvcc 编译 CUDA 程序

简单 CUDA 程序的基本框架

单源文件 CUDA 程序基本框架

自定义设备函数

CUDA 程序的错误检测

检测 CUDA 运行时错误的宏函数

CUDA-MEMCHECK 检查内存错误

获得 GPU 加速的关键

CUDA 事件计时

nvprof 查看程序性能

影响 GPU 加速的关键因素

CUDA 的数学函数库

CUDA 的内存组织

CUDA 中不同类型的内存

全局内存

常量内存

纹理内存

寄存器

局部内存

共享内存

L1 和 L2 缓存

SM 及其占有率

CUDA 运行时 API 函数查询设备

全局内存的合理使用

全局内存的合并与非合并访问

矩阵转置

共享内存的合理使用

数组归约

矩阵转置

共享内存的 bank 冲突

原子函数的合理使用

完全在 GPU 中进行归约

原子函数

邻居列表

线程束基本函数与协作组

单指令-多线程模式

线程束内的线程同步函数

更多线程束内的基本函数

协作组

数组归约程序的进一步优化

CUDA 流

在默认流中重叠主机和设备计算

非默认 cuda 流重叠多个核函数

非默认 cuda 流重叠核函数与数据传递

分子动力学模型

cuda 标准库的使用

Thrust

cuBLAS

cuSolver

cuRAND

1
2
3
>> g++ hello.cpp -o ./bin/hello.exe
>> ./bin/hello
msvc: hello world!

使用核函数的 CUDA 程序