NVIDA CUDA-DirverAPI入门

NVIDA公司发布了CUDA，CUDA是建立在NVIDA的CPUs上的一个通用并行计算平台和编程模型，基于CUDA可以利用GPUs的并行计算引擎来更加高效地解决比较复杂的计算难题。近年来，GPU最成功的应用领域就是深度学习领域，基于GPU的并行计算已经成为训练深度学习模型的标配了。

GPU并不是一个独立运行的计算平台，需要与CPU协同工作，可以看作CPU的协处理器，因此当我们在说GPU在并行计算时，其实是基于CPU+GPU的异构计算架构。在异构计算架构中，GPU和CPU通过PCle总线连接在一起来协同工作，CPU所在位置称为主机端host，而GPU所在位置称为设备端device。

CUDA是NVIDA公司所开发的GPU编程模型，它提供了GPU编程的简易接口，基于CUDA编程可以构建基于GPU计算的应用程序。

0x01 CUDA编程基础

（一）CUDA认知

（二）显卡、显卡驱动、NVCC、CUDA driver、CUDAtoolkit、cudnn是什么？

（三）runtime和driver API的区别

0x02 环境配置

0x03 cuInit——驱动初始化

0x04 CUcontext——上下文管理

（一）简介

（二）使用代码

0x05 Memory

（一）Pinned Memory与Pageable Memory

（一）编写核函数必须遵循CUDA规范，那么有哪些规范？

（二）函数声明中，__ global __ 、__ device __ 、__ host __三者区别是什么？

（三）核函数内部怎么使用线程？为什么要分成三个层次？

0x0B Thrust and Error

0x01 CUDA编程基础

我觉得参考价值较大的文章：

CUDA编程入门极简教程 - 知乎

显卡，显卡驱动,nvcc, cuda driver,cudatoolkit,cudnn到底是什么？ - 知乎

（一）CUDA认知

CUDA是一个异构模型，需要CPU和GPU协同工作。
在CUDA中，host指代CPU及其内存，device指代GPU及其内存。
CUDA程序中既包含host程序又包含device程序，他们分别在CPU和GPU上运行，同时CPU与GPU之间还可以进行通讯。
CUDA Driver是与GPU沟通的驱动级别底层API，对应于cuda.h和libcuda.so文件
典型的CUDA程序的执行流程如下：
- 分配host内存，并进行数据初始化。
- 分配device内存，并从host将数据拷贝到device上。
- 调用CUDA的核函数在device上完成指定的运算。
- 将device上的运算结果拷贝到host上。
- 释放device和host上分配的内存。

在上面的过程中，调用了CUDA的核函数来执行并行计算，kernel是CUDA中一个重要的概念，kernel是在device上线程中并执行的函数，核函数用__ global __符号声明，在调用中需要用<<grid,block>>来指定kernel要执行的线程数量。在CUDA中，每一个线程都要执行核函数，并且每个线程会分配一个唯一的线程号threadID，这个ID值可以通过核函数的内置变量threadIdx来获得。

由于GPU实际上是异构模型，所以需要区分host和device上的代码，在CUDA中是通过函数类型限定词来区分开host和device上的函数，主要函数类型限定词如下：

__ global __ ：在device上执行，从host中调用（一些特定的GPU也可以从device中调用），返回类型必须是void，不支持可变参数，不能成为类成员函数。注意用 __ global __定义的kernel是异步的，这意味着host不会等待kernel执行完就执行下一步。

__ device __ ：在device上执行，但仅可以从device中调用，不可以和 __ global __ 同时用。

__ host __ ：在host上执行，仅可以从host上调用，一般省略不写，不可以和 __ global __同时用，但可和 _ _ device _ _ ，此时函数会在device和host都编译。

（二）显卡、显卡驱动、NVCC、CUDA driver、CUDAtoolkit、cudnn是什么？

显卡：简单理解这个就是我们前面说的GPU，尤其指NVIDIA公司生产的GPU系列，因为后面介绍的cuda,cudnn都是NVIDIA公司针对自身的GPU独家设计的。
显卡驱动：驱动软件，类比声卡驱动，摄像头驱动。很明显就是字面意思，通常指NVIDIA Driver，其实它就是一个驱动软件，而前面的显卡就是硬件。
GPU架构：硬件设计方式。现在有Tesla、Fermi、Kepler、Maxwell、Pascal等。
CUDA：其中一种理解是它是一种编程语言（像c++,python等，只不过它是专门用来操控GPU的）
cudnn：这个其实就是一个专门为深度学习计算设计的软件库，里面提供了很多专门的计算函数，它里面包含了许多库，例如：cudart, cublas等。
NVCC：nvcc其实就是CUDA的编译器,可以从CUDA Toolkit的/bin目录中获取,类似于gcc就是c语言的编译器。由于程序是要经过编译器编程成可执行的二进制文件，而cuda程序有两种代码，一种是运行在cpu上的host代码，一种是运行在gpu上的device代码，所以nvcc编译器要保证两部分代码能够编译成二进制文件在不同的机器上执行。
nvidia -smi：nvidia-smi全程是NVIDIA System Management Interface ，它是一个基于前面介绍过的NVIDIA Management Library(NVML)构建的命令行实用工具，旨在帮助管理和监控NVIDIA GPU设备。

那么NVCC、nvidia -smi的区别在哪：

我们在windows上输入nvcc --version以及nvidia-smi得到的CUDA版本信息是不一样的。

>nvcc --version
nvcc: NVIDIA (R) Cuda compiler driver
Copyright (c) 2005-2021 NVIDIA Corporation
Built on Fri_Dec_17_18:28:54_Pacific_Standard_Time_2021
Cuda compilation tools, release 11.6, V11.6.55
Build cuda_11.6.r11.6/compiler.30794723_0

+-----------------------------------------------------------------------------+
| NVIDIA-SMI 511.23       Driver Version: 511.23       CUDA Version: 11.6     |
|-------------------------------+----------------------+----------------------+
| 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  NVIDIA GeForce ... WDDM  | 00000000:01:00.0 Off |                  N/A |
| N/A   56C    P8    N/A /  N/A |      0MiB /  4096MiB |      1%      Default |
|                               |                      |                  N/A |
+-------------------------------+----------------------+----------------------+

+-----------------------------------------------------------------------------+
| Processes:                                                                  |
|  GPU   GI   CI        PID   Type   Process name                  GPU Memory |
|        ID   ID                                                   Usage      |
|=============================================================================|
|  No running processes found                                                 |
+-----------------------------------------------------------------------------+

CUDA有两个主要的API：runtime(运行时) API和driver API。这两个API都有对应的CUDA版本（runtime API和driver API），在网上的某个解释是这样的：

来源于：Different CUDA versions shown by nvcc and NVIDIA-smi - Stack Overflow

用于支持driver API的必要文件(如libcuda.so)是由GPU driver installer安装的。nvidia-smi就属于这一类API。
用于支持runtime API的必要文件(如libcudart.so以及nvcc)是由CUDA Toolkit installer安装的。（CUDA Toolkit Installer有时可能会集成了GPU driver Installer）。nvcc是与CUDA Toolkit一起安装的CUDA compiler-driver tool，它只知道它自身构建时的CUDA runtime版本。它不知道安装了什么版本的GPU driver，甚至不知道是否安装了GPU driver。

综上，如果driver API和runtime API的CUDA版本不一致可能是因为你使用的是单独的GPU driver installer，而不是CUDA Toolkit installer里的GPU driver installer。

（三）runtime和driver API的区别

runtime和driverAPI在很多情况下非常相似，用起来其实都是差不多的，但是这两个API不可以混用，两者是互斥的，也就是说在开发时我们只能选择其中一种API。

简单的二者的区别：runtime是更高级别的封装，开发人员用起来更方便，而driver API更接近底层，速度可能会更快。

CUDA Driver与CUDA Runtime相比更偏底层，就意味着Driver API有着更灵活的控制，也伴随着更复杂的编程。
因此CUDA driver需要做显式的初始化cuInit(0)，否则其他API都会返回CUDA_ERROR_NOT_INITIALIZED

0x02 环境配置

这里网上有很多教程，那么就不过多描述：

0x03 cuInit——驱动初始化

cuInit的意义是：初始化驱动API，如果不执行，则所有API都将返回错误，全局执行一次即可。

cuInit没有对应的cuDestroy，不需要释放，程序销毁自动释放。

cuInit的参数以及函数说明可以在NVIDA官网中查看：

CUDA Driver API :: CUDA Toolkit Documentation

说明：

 cuInit(int flags), 这里的flags目前必须给0;
 对于cuda的所有函数，必须先调用cuInit，否则其他API都会返回CUDA_ERROR_NOT_INITIALIZED

那么我们就先实现使用CUDA驱动头文件来进行初始化，获取我们的CUDA版本及其获取当前设备信息吧：

获取驱动版本管理：CUDA 驱动程序 API ：： CUDA Toolkit Documentation (nvidia.com)

获取设备信息管理：CUDA 驱动程序 API ：： CUDA Toolkit Documentation (nvidia.com)

// CUDA驱动头文件cuda.h
#include <cuda.h>
#include <stdio.h> 
#include <string.h>
int main(){

    /* 
    cuInit(int flags), 这里的flags目前必须给0;
    对于cuda的所有函数，必须先调用cuInit，否则其他API都会返CUDA_ERROR_NOT_INITIALIZED
     */
    CUresult code=cuInit(0);  //CUresult 类型：用于接收一些可能的错误代码
    if(code != CUresult::CUDA_SUCCESS){
        const char* err_message = nullptr;
        cuGetErrorString(code, &err_message);    // 获取错误代码的字符串描述
        // cuGetErrorName (code, &err_message);  // 也可以直接获取错误代码的字符串
        printf("Initialize failed. code = %d, message = %s\n", code, err_message);
        return -1;
    }
    /* 
    测试获取当前cuda驱动的版本
    显卡、CUDA、CUDA Toolkit

    1. 显卡驱动版本，比如：Driver Version: 460.84
    2. CUDA驱动版本：比如：CUDA Version: 11.2
    3. CUDA Toolkit版本：比如自行下载时选择的10.2、11.2等；这与前两个不是一回事, CUDA Toolkit的每个版本都需要最低版本的CUDA驱动程序
        
    三者版本之间有依赖关系, 可参照https://docs.nvidia.com/cuda/cuda-toolkit-release-notes/index.html
    nvidia-smi显示的是显卡驱动版本和此驱动最高支持的CUDA驱动版本
        
     */
    int driver_version = 0;
    code = cuDriverGetVersion(&driver_version);  // 获取驱动版本
    printf("CUDA Driver version is %d\n", driver_version); // 若driver_version为11020指的是11.2

    // 测试获取当前设备信息
    char device_name[100]; // char 数组
    CUdevice device = 0;
    code = cuDeviceGetName(device_name, sizeof(device_name), device);  // 获取设备名称、型号如：Tesla V100-SXM2-32GB // 数组名device_name当作指针
    printf("Device %d name is %s\n", device, device_name);
    return 0;
}

那么NVIDA中如何对没有初始化的错误应该如何检查，NVIDA官网提供的查错函数：

// 使用有参宏定义检查cuda driver是否被正常初始化, 并定位程序出错的文件名、行数和错误信息
// 宏定义中带do...while循环可保证程序的正确性
#define checkDriver(op)    \
    do{                    \
        auto code = (op);  \
        if(code != CUresult::CUDA_SUCCESS){     \
            const char* err_name = nullptr;     \
            const char* err_message = nullptr;  \
            cuGetErrorName(code, &err_name);    \
            cuGetErrorString(code, &err_message);   \
            printf("%s:%d  %s failed. \n  code = %s, message = %s\n", __FILE__, __LINE__, #op, err_name, err_message);   \
            return -1;   \
        }                \
    }while(0)

那么我们首先就来使用一下，相对于是一个helloword吧：

// CUDA驱动头文件cuda.h
#include <cuda.h>

#include <stdio.h>
#include <string.h>

// 使用有参宏定义检查cuda driver是否被正常初始化, 并定位程序出错的文件名、行数和错误信息
// 宏定义中带do...while循环可保证程序的正确性
#define checkDriver(op)    \
    do{                    \
        auto code = (op);  \
        if(code != CUresult::CUDA_SUCCESS){     \
            const char* err_name = nullptr;     \
            const char* err_message = nullptr;  \
            cuGetErrorName(code, &err_name);    \
            cuGetErrorString(code, &err_message);   \
            printf("%s:%d  %s failed. \n  code = %s, message = %s\n", __FILE__, __LINE__, #op, err_name, err_message);   \
            return -1;   \
        }                \
    }while(0)

int main(){

    //检查cuda driver的初始化。虽然不初始化或错误初始化某些API不会报错（不信你试试），但安全起见调用任何API前务必检查cuda driver初始化
    cuInit(2); // 正确的初始化应该给flag = 0
    checkDriver(cuInit(0));

    // 测试获取当前cuda驱动的版本
    int driver_version = 0;
    checkDriver(cuDriverGetVersion(&driver_version));
    printf("Driver version is %d\n", driver_version);

    // 测试获取当前设备信息
    char device_name[100];
    CUdevice device = 0;
    checkDriver(cuDeviceGetName(device_name, sizeof(device_name), device));
    printf("Device %d name is %s\n", device, device_name);
    return 0;
}

其实官方给的查错代码还是有很多值得改进的地方的，首先它使用宏定义，每当我需要再加一行的时候，就需要再加一个“\”，这样很麻烦，而且可以看到do...while(0)的使用与顺序执行...的效果是一样的，但是前者可以保证程序的正确性，下面提供了一种改进的方法：

// 很明显，这种代码封装方式，更加的便于使用
//宏定义 #define <宏名>（<参数表>） <宏体>
#define checkDriver(op)  __check_cuda_driver((op), #op, __FILE__, __LINE__)

bool __check_cuda_driver(CUresult code, const char* op, const char* file, int line){

    if(code != CUresult::CUDA_SUCCESS){    
        const char* err_name = nullptr;    
        const char* err_message = nullptr;  
        cuGetErrorName(code, &err_name);    
        cuGetErrorString(code, &err_message);   
        printf("%s:%d  %s failed. \n  code = %s, message = %s\n", file, line, op, err_name, err_message);   
        return false;
    }
    return true;
}

那么我们的查错程序就可以这么修改：

// CUDA驱动头文件cuda.h
#include <cuda.h>
#include <stdio.h>
#include <string.h>

// 很明显，这种代码封装方式，更加的便于使用
//宏定义 #define <宏名>（<参数表>） <宏体>
#define checkDriver(op)  __check_cuda_driver((op), #op, __FILE__, __LINE__)

bool __check_cuda_driver(CUresult code, const char* op, const char* file, int line){

    if(code != CUresult::CUDA_SUCCESS){    
        const char* err_name = nullptr;    
        const char* err_message = nullptr;  
        cuGetErrorName(code, &err_name);    
        cuGetErrorString(code, &err_message);   
        printf("%s:%d  %s failed. \n  code = %s, message = %s\n", file, line, op, err_name, err_message);   
        return false;
    }
    return true;
}

int main(){

    // 检查cuda driver的初始化
    // 实际调用的是__check_cuda_driver这个函数
    checkDriver(cuInit(0));

    // 测试获取当前cuda驱动的版本
    int driver_version = 0;
    if(!checkDriver(cuDriverGetVersion(&driver_version))){
        return -1;
    }
    printf("Driver version is %d\n", driver_version);

    // 测试获取当前设备信息
    char device_name[100];
    CUdevice device = 0;
    checkDriver(cuDeviceGetName(device_name, sizeof(device_name), device));
    printf("Device %d name is %s\n", device, device_name);
    return 0;
}

0x04 CUcontext——上下文管理

（一）简介

参考：交互式 GPU 编程 - 第 3 部分 - CUDA 上下文恶作剧 (dragan.rocks)

context是一种上下文，关联对GPU的所有操作。设备与特定进程相关连的所有状态。比如，你写的一段kernel code对GPU的使用会造成不同状态（内存映射、分配、加载的code），Context则保存着所有的管理数据来控制和使用设备。
context与一块显卡关联，一个显卡可以被多个context关联。GPU的context相对于CPU的program，一块GPU上可以有多个contexts，但是它们之间是相互隔离的，建议一块设备就一个context。
每个线程都有一个栈结构储存context，栈顶是当前使用的context，对应有push、pop函数操作context的栈，所有api都以当前context为操作目标。
它的作用是什么呢，我们可以想象到我们每次进行任何操作都需要传递一个device决定送到哪个设备执行，很麻烦，可以看看下面这种没有context的代码：

那么再看看有context的代码：

这样可以说是更安全且更方便了吧
context只是为了方便控制device的一种手段而提出来的，栈的存在是为了方便控制多个设备。
由于高频操作，是一个线程基本固定访问一个显卡不变，且只使用一个context，很少会用到多context。
CreateContext、PushCurrent、PopCurrent这种多context的管理就显得很麻烦了，还可以再简单，因此推出了cuDevicePrimaryCtxRetain，为设备关联主context，分配、释放、设置、栈都不用你来控制管理了：

看看有context的代码：

再看看不用管理context栈的代码：（runtimeAPI自动使用cuDevicePrimaryCtxRetain）

primaryContext：一个显卡对应一个primary context，你只需要提供一个设备id，就可以把你的context设置好。
不同线程，只要设备id一样，primary context就一样，context是线程安全的。

上下文可以做什么事情呢：

持有分配的内存列表
持有加载进设备的kernel code
CPU与GPU之间的unified memory
...

如何管理上下文？

在cuda driver同样需要显示管理上下文

开始时cuCtxCreate()创建上下文，结束时cuCtxDestroy销毁上下文。像文件管理一样须手动开关。

用cuDevicePrimaryCtxRetain()创建上下文更好！

cuCtxGetCurrent()获取当前上下文

可以使用堆栈管理多个上下文cuCtxPushCurrent()压入，cuCtxPopCurrent()推出

对ctxA使用cuCtxPushCurrent()和cuCtxCreate()都相当于将ctxA放到栈顶（让它成为current context）

cuda runtime可以自动创建，是基于cuDevicePrimaryCtxRetain()创建的。

（二）使用代码

// CUDA驱动头文件cuda.h
#include <cuda.h>   // include <> 和 "" 的区别    
#include <stdio.h>  // include <> : 标准库文件 
#include <string.h> // include "" : 自定义文件  

#define checkDriver(op)  __check_cuda_driver((op), #op, __FILE__, __LINE__)

bool __check_cuda_driver(CUresult code, const char* op, const char* file, int line){
    if(code != CUresult::CUDA_SUCCESS){    // 如果 成功获取CUDA情况下的返回值 与我们给定的值(0)不相等， 即条件成立， 返回值为flase
        const char* err_name = nullptr;    // 定义了一个字符串常量的空指针
        const char* err_message = nullptr;  
        cuGetErrorName(code, &err_name);    
        cuGetErrorString(code, &err_message);   
        printf("%s:%d  %s failed. \n  code = %s, message = %s\n", file, line, op, err_name, err_message); //打印错误信息
        return false;
    }
    return true;
}

int main(){

    // 检查cuda driver的初始化
    checkDriver(cuInit(0));

    // 为设备创建上下文
    CUcontext ctxA = nullptr;                                   // CUcontext 其实是 struct CUctx_st*（是一个指向结构体CUctx_st的指针）
    CUcontext ctxB = nullptr;
    CUdevice device = 0;
    checkDriver(cuCtxCreate(&ctxA, CU_CTX_SCHED_AUTO, device)); // 这一步相当于告知要某一块设备上的某块地方创建 ctxA 管理数据。输入参数 参考 https://www.cs.cmu.edu/afs/cs/academic/class/15668-s11/www/cuda-doc/html/group__CUDA__CTX_g65dc0012348bc84810e2103a40d8e2cf.html
    checkDriver(cuCtxCreate(&ctxB, CU_CTX_SCHED_AUTO, device)); // 参考 1.ctx-stack.jpg
    printf("ctxA = %p\n", ctxA);
    printf("ctxB = %p\n", ctxB);
    /* 
        contexts 栈：
            ctxB -- top <--- current_context
            ctxA 
            ...
     */

    // 获取当前上下文信息
    CUcontext current_context = nullptr;
    checkDriver(cuCtxGetCurrent(&current_context));             // 这个时候current_context 就是上面创建的context
    printf("current_context = %p\n", current_context);

    // 可以使用上下文堆栈对设备管理多个上下文
    // 压入当前context
    checkDriver(cuCtxPushCurrent(ctxA));                        // 将这个 ctxA 压入CPU调用的thread上。专门用一个thread以栈的方式来管理多个contexts的切换
    checkDriver(cuCtxGetCurrent(&current_context));             // 获取current_context (即栈顶的context)
    printf("after pushing, current_context = %p\n", current_context);
    /* 
        contexts 栈：
            ctxA -- top <--- current_context
            ctxB
            ...
    */
    
    // 弹出当前context
    CUcontext popped_ctx = nullptr;
    checkDriver(cuCtxPopCurrent(&popped_ctx));                   // 将当前的context pop掉，并用popped_ctx承接它pop出来的context
    checkDriver(cuCtxGetCurrent(&current_context));              // 获取current_context(栈顶的)
    printf("after poping, popped_ctx = %p\n", popped_ctx);       // 弹出的是ctxA
    printf("after poping, current_context = %p\n", current_context); // current_context是ctxB

    checkDriver(cuCtxDestroy(ctxA));
    checkDriver(cuCtxDestroy(ctxB));

    // 更推荐使用cuDevicePrimaryCtxRetain获取与设备关联的context
    // 注意这个重点，以后的runtime也是基于此, 自动为设备只关联一个context
    checkDriver(cuDevicePrimaryCtxRetain(&ctxA, device));       // 在 device 上指定一个新地址对ctxA进行管理
    printf("ctxA = %p\n", ctxA);
    checkDriver(cuDevicePrimaryCtxRelease(device));
    return 0;
}

0x05 Memory

首先要学会内存的各个单元，这样才可以使用cuda进行高效的内存分配。

CUDA的内存模型：

首先可以看看GPU的电路板：

那么GPU与电脑主板的关系：

（一）Pinned Memory与Pageable Memory

对于整个Host Memory内存条而言，操作系统区分为两个大类（逻辑区分，物理上是同一个东西）：

Pageable Memory：可分页内存
Page Lock Memory（Pinned Memory）：页锁定内存

怎么理解呢：

Pageable Memory可以理解为是VIP房间，锁定给你一个人用。而Page Lock Memory是普通房间，当酒店的房间不够用的时候，会选择性把你的房间腾出来给其他人用，这就可以容纳更多人了。这其实有点像是虚拟内存的东西，造成了房间很多的这个假象，代价是性能降低。

那么再做一下总结：

pinned memory具有锁定特性，是稳定不会被交换的（这很重要，相当于每次去这个房间都一定能找到你）。
pageable memory没有锁定特性，对于第三方设备（比如GPU），去访问时，因为无法感知内存是否被交换，可能得不到正确的数据（每次去房间找，说不准你的房间被人交换了）。
pageable memory的性能比pinned memory差，很可能降低你程序的优先级然后把内存交换给别人用。
pageable memory策略能使用内存假象，实际8GB但是可以使用15GB，提高程序运行数量（不是速度）。
pinned memory太多，会导致操作系统整体性能降低（程序运行数量减少），8GB就只能用8GB。注意不是你的应用程序性能降低。
GPU可以直接访问pinned memory而不能访问pageable memory。

那么当数据传到GPU时，路径是这样的：

当我们使用Pageable Memory时，我们还需要经过Pinned Memory才可以把数据传到DRAM，这当我们数据规模大时，大打折扣；而使用Pinned Memory时，我们是直接传输到DRAM。也就是说当我们使用显卡访问Pinnned Memory轨迹时，是这样的：

通过PICE接口，到主板，再到内存条。

那么我们再详细地概括一下总结：

1.GPU可以直接访问pinned memory，称之为（DMA Direct Memory Access），以前接触单片机的时候觉得DMA只是一个越过CPU访问内存的一个小东西，没想到在这里直接上升为GPU。

2.对于GPU访问而言，距离计算单元越近，效率越高，所以PinnedMemory<GlobalMemory<SharedMemory

3.代码中，由new、malloc分配的，是pageable memory，由cudaMallocHost分配的是PinnedMemory，由cudaMalloc分配的是GlobalMemory。

4.尽量多用PinnedMemory储存host数据，或者显式处理Host到Device时，用PinnedMemory做缓存，都是提高性能的关键。

（二）内存分配

统一内存（Unified addressing）：

参考博客：初步介绍CUDA中的统一内存_扫地的小何尚的博客-CSDN博客_cuda统一内存

可以非常轻松地分配合访问可由系统中任何处理器、CPU或GPU上运行代码使用的数据。这种硬件/软件技术允许应用程序分配可以从CPU或GPU上运行的代码读取或写入的数据。

分配线性内存cuMemAlloc():

线性内存：线性内存被组织在单个连续的地址空间中，可以直接以及线性地访问这些内存位置。内存分配空间以字节为大小，并返回所分配的内存地址。

分配主机锁页内存cuMemAllocHost():

锁页内存：页面不允许被调入调出的叫锁页内存，反之叫可分页内存。

好处：快。

a. 设备可以直接访问内存，与可分页内存相比，它的读写带宽要高得多
b. 驱动程序会跟踪使用cuMemAllocHost()分配的虚拟内存范围，并自动加速对cuMemcpy（）等函数的调用。
使用注意：分配过多锁业内存会减少系统可用于分页的内存量，可能会降低系统性能。因此，在主机和设备之间为数据交换分配临时区域时，最好少用此功能。

这里是从主机分配内存，因此不是输入device prt的地址，而是一个主机的二级地址。

使用：

内存的初始化cuMemsetD32(CUdeviceptr dstDevice, unsigned int ui, size_t N), 将N个32位值的内存范围设置为指定的值ui。

内存的释放cuMemFreeHost(): 有借有还再借不难。

选择右边的内存分配模型：

参考博客：用于数据传输的页面锁定主机内存 - 雷毛的日志 (leimao.github.io)

主机和 CUDA 设备之间的数据传输需要主机上的页面锁定内存。如果数据不在主机上的可分页内存中，则在将数据从主机传输到设备期间，数据将从可分页主机内存隐式传输到临时页面锁定的主机内存，然后将数据从页面锁定的主机内存传输到设备内存。

这会带来额外的数据传输开销，并且在某些情况下可能会显著影响计算机程序的性能。因此，问题是为什么主机和CUDA设备之间的数据传输需要在主机上使用页面锁定内存。我认为主要动机是保证数据传输效率。

由于主机和CUDA设备之间的数据传输无论如何都必须使用页面锁定的内存，因此为了优化数据传输，我们可以直接分配页面锁定内存来存储数据，而不是分配可分页内存来存储数据。

需要注意的是，手动使用大量页面锁定内存可能会导致操作系统性能问题。手动使用页面锁定内存意味着用户负责分配和释放页面锁定内存（内存泄漏）。如果用户未及时释放页面锁定内存，则物理内存上可用于新数据和新应用程序的可用内存将减少。因此，操作系统可能会变得不稳定。如果页面锁定内存仅在可分页内存和 CUDA 内存之间的数据传输期间临时使用，则它将始终及时释放，但代价是在可分页内存和页面锁定内存之间增加数据传输开销。

因此，我们被鼓励使用页面锁定的内存，但我们不应该滥用它。通常，如果我们知道数据将在主机和CUDA设备之间多次传输，则最好将数据放在页面锁定的内存中以避免不必要的开销。

（三）代码

// CUDA驱动头文件cuda.h
#include <cuda.h>
#include <stdio.h>
#include <string.h>

#define checkDriver(op)  __check_cuda_driver((op), #op, __FILE__, __LINE__)

bool __check_cuda_driver(CUresult code, const char* op, const char* file, int line){

    if(code != CUresult::CUDA_SUCCESS){    
        const char* err_name = nullptr;    
        const char* err_message = nullptr;  
        cuGetErrorName(code, &err_name);    
        cuGetErrorString(code, &err_message);   
        printf("%s:%d  %s failed. \n  code = %s, message = %s\n", file, line, op, err_name, err_message);   
        return false;
    }
    return true;
}

int main(){

    // 检查cuda driver的初始化
    checkDriver(cuInit(0));

    // 创建上下文
    CUcontext context = nullptr;
    CUdevice device = 0;
    checkDriver(cuCtxCreate(&context, CU_CTX_SCHED_AUTO, device));
    printf("context = %p\n", context);

    // 输入device prt向设备要一个100 byte的线性内存，并返回地址
    CUdeviceptr device_memory_pointer = 0;
    checkDriver(cuMemAlloc(&device_memory_pointer, 100)); // 注意这是指向device的pointer, 
    printf("device_memory_pointer = %p\n", device_memory_pointer);

    // 输入二级指针向host要一个100 byte的锁页内存，专供设备访问。
    float* host_page_locked_memory = nullptr;
    checkDriver(cuMemAllocHost((void**)&host_page_locked_memory, 100));
    printf("host_page_locked_memory = %p\n", host_page_locked_memory);

    // 向page-locked memory 里放数据（仍在CPU上），可以让GPU可快速读取
    host_page_locked_memory[0] = 123;
    printf("host_page_locked_memory[0] = %f\n", host_page_locked_memory[0]);
    /* 
        记住这一点
        host page locked memory 声明的时候为float*型，可以直接转换为device ptr，这才可以送给cuda核函数（利用DMA(Direct Memory Access)技术）
        初始化内存的值: cuMemsetD32 ( CUdeviceptr dstDevice, unsigned int  ui, size_t N )
        初始化值必须是无符号整型，因此需要将new_value进行数据转换：
        但不能直接写为:(int)value，必须写为*(int*)&new_value, 我们来分解一下这条语句的作用：
        1. &new_value获取float new_value的地址
        2.(int*)将地址从float * 转换为int*以避免64位架构上的精度损失
        3.*(int*)取消引用地址，最后获取引用的int值
     */
    
    float new_value = 555;
    checkDriver(cuMemsetD32((CUdeviceptr)host_page_locked_memory, *(int*)&new_value, 1)); 
    printf("host_page_locked_memory[0] = %f\n", host_page_locked_memory[0]);

    // 释放内存
  checkDriver(cuMemFreeHost(host_page_locked_memory));
    return 0;
}

0x06 stream - 流

CUDA流表示一个GPU操作队列，该队列中的操作将以添加到流中的先后顺序而一次执行。可以将一个流看做是GPU上的一个任务，不同的任务可以并行执行。使用CUDA流，首先要选择一个支持设备重叠功能的设备，支持设备重叠功能的GPU能够在执行一个CUDA核函数的同时，还能在主机和设备之间执行复制数据操作。

支持重叠功能的设备的这一特性很重要，可以在一定程度上提升GPU程序的执行效率。一般情况下，CPU内存远大于GPU内存，对于数据量比较大的情况，不可能把CPU缓冲区中的数据一次性传输给GPU，需要分块传输，如果能够在分块传输的同时，GPU也在执行核函数运算，这样的异步操作，就用到设备的重叠功能，能够提高运算性能。说白了，其实就是这个设备具有多个核，可以同时进行工作。

流是一种基于context之上的任务管道抽象，一个context可以创建n个流。
流是异步控制的主要方式。
nullptr表示默认流，每个线程都有自己的默认流。
某一方发出指令后，他就可以做任何事情，无需等待指令执行完毕，指令发出的耗时也是极短的，这是一个异步操作，执行的代码加入流的队列后，立即返回，不耽误时间。
接收指令的一方，根据流的队列，顺序执行；发出指令的一方具有选择性，在需要的时候等待所有的执行结果。
新建一个流，就是新建一个接收指令的一方，可以新建很多个流。
通过cudaEvent可以选择性等待任务队列中的部分任务是否就绪。

需要注意的是：

要十分注意，指令发出后，流队列中储存的是指令参数，不能加入队列后立即释放参数指针，这会导致流队列执行该指令时指针失效而出错。
应当在十分肯定流已经不需要这个指针后，才进行修改或者释放，否则会有非预期结果出现。

总结来说：

stream是一个流句柄，可以当做是一个队列
- cuda执行器从stream中一条条的读取并执行指令。
- 例如cudaMemcpyAsync函数等同于向stream这个队列中加入一个cudaMemcpy指令并排队。
- 使用到了stream的函数，便立即向stream中加入指令后立即返回，并不会等待指令执行结束。
- 通过cudaStreamSynchronize函数，等待stream中所有指令执行完毕，也就是队列为空。
当使用stream时，需要注意
- 由于异步函数会立即返回，因此传递进入的参数要考虑其生命周期，应确认函数调用结束后再做释放。
还可以向stream中加入Event，用以监控是否到达了某个检查点
- cudaEventCreate，创建事件
- cudaEventRecord，记录事件，即在stream中加入某个事件，当队列执行到该事件后，修改其状态
- cudaEventQuery，查询事件当前状态
- cudaEventElapsedTime，计算两个事件之间经历的时间间隔，若要统计某些核函数执行时间，请使用这个函数，能够得到最准确的统计
- cudaEventSynchronize，同步某个事件，等待事件到达
- cudaStreamWaitEvent，等待流中的某个事件
默认流，对于cudaMemcpy等同步函数，其等价于执行了
- cudaMemcpyAsync(... 默认流) 加入队列
- cudaStreamSynchronize(默认流) 等待执行完成
- 默认流与当前设备上下文类似，是与当前设备进行的关联
- 因此，如果大量使用默认流，会导致性能低下

如何使用：

// CUDA运行时头文件
#include <cuda_runtime.h>
#include <stdio.h>
#include <string.h>

#define checkRuntime(op)  __check_cuda_runtime((op), #op, __FILE__, __LINE__)

bool __check_cuda_runtime(cudaError_t code, const char* op, const char* file, int line){
    if(code != cudaSuccess){    
        const char* err_name = cudaGetErrorName(code);    
        const char* err_message = cudaGetErrorString(code);  
        printf("runtime error %s:%d  %s failed. \n  code = %s, message = %s\n", file, line, op, err_name, err_message);   
        return false;
    }
    return true;
}

int main(){

    int device_id = 0;
    checkRuntime(cudaSetDevice(device_id));

    cudaStream_t stream = nullptr;
    checkRuntime(cudaStreamCreate(&stream));

    // 在GPU上开辟空间
    float* memory_device = nullptr;
    checkRuntime(cudaMalloc(&memory_device, 100 * sizeof(float)));

    // 在CPU上开辟空间并且放数据进去，将数据复制到GPU
    float* memory_host = new float[100];
    memory_host[2] = 520.25;
    checkRuntime(cudaMemcpyAsync(memory_device, memory_host, sizeof(float) * 100, cudaMemcpyHostToDevice, stream)); // 异步复制操作，主线程不需要等待复制结束才继续

    // 在CPU上开辟pin memory,并将GPU上的数据复制回来 
    float* memory_page_locked = nullptr;
    checkRuntime(cudaMallocHost(&memory_page_locked, 100 * sizeof(float)));
    checkRuntime(cudaMemcpyAsync(memory_page_locked, memory_device, sizeof(float) * 100, cudaMemcpyDeviceToHost, stream)); // 异步复制操作，主线程不需要等待复制结束才继续
    checkRuntime(cudaStreamSynchronize(stream));
    
    printf("%f\n", memory_page_locked[2]);
    
    // 释放内存
    checkRuntime(cudaFreeHost(memory_page_locked));
    checkRuntime(cudaFree(memory_device));
    checkRuntime(cudaStreamDestroy(stream));
    delete [] memory_host;
    return 0;
}

0x07 核函数

CUDA核函数是指GPU端运行的代码，核函数内部主要干了什么？简而言之，就是规定GPU的各个线程访问哪个数据并执行什么计算。

通过xxx.cu创建一个cudac程序文件，并把cu交给nvcc编译，才能识别cuda语法。

参考博客：

1.2.CUDA核函数 - 知乎

（一）编写核函数必须遵循CUDA规范，那么有哪些规范？

必须写在*.cu文件中。
必须以__ global __限定符声明定义（只有他才可以执行核函数）。
返回类型必须是void。
不支持可变数量的参数。调用核函数是传值的，不能传引用，可以传递类、结构体等，核函数可以是模板。
核函数内部只能访问设备内存。
核函数内部不能使用静态变量。

（二）函数声明中， global 、 device 、 host 三者区别是什么？

1） __ global __ 修饰的函数是核函数，在设备端执行，可以从主机端调用，也可以在sm3以上的设备端调用（比如动态并行）；只有__ global __ 修饰的函数才可以用<<<>>>的方式调用。 2） __ device __ 修饰的函数是设备函数，在设备端执行，只能从设备端调用； 3） __ host __ 修饰的函数是主机函数，在主机端执行，只能从主机端调用；

4） __ device __ 和 __ host __ 可以一起使用，来表示该函数可以同时在主机端和设备端执行；

5）__ shared __ 表示变量为共享变量。

6） nvcc编译选项中添加-dc（相当于--relocatable-device-code=true --compile）时，global函数可以调用其它文件中的device函数，否则只能调用同文件中的device函数。

（三）核函数内部怎么使用线程？为什么要分成三个层次？

CUDA从逻辑上将GPU线程分成了三个层次——线程格（grid）、线程块（block）和线程（thread）。

每个核函数对应一个线程格，一个线程格中有一个或多个线程块，一个线程块中有一个或多个线程。

关于GPU的硬件描述可以看这篇：

2.1.GPU硬件架构 - 知乎

CUDA核函数中文说明要将线程分为三个层次，其实是与GPU的硬件组成相关联的。在GPU硬件中本身就存在三个层次——核心、流多处理器、设备，这是一种类似于计算机集群的层次结构，而我们编写的核函数正是运行在这种层次结构上，所以核函数必须支持这三个层次，否则任务无法顺利分解，也就无法从高层次向低层次传递。

实际设计中，CUDA将这种对应关系规定为：

Grid分配到Device上运行。
Block分配到SM上运行。SM可以理解为一个计算机集群。
Thread分配到Croe上运行。Core为计算核心，负责浮点数和整数计算。

（四）cu文件引入的新的符号和语法强调

__global__标记，核函数标记
- 调用方必须是host
- 返回值必须是void
- 例如：__global__ void kernel(const float* pdata, int ndata)
- 必须以kernel<<<gridDim, blockDim, bytesSharedMemorySize, stream>>>(pdata, ndata)的方式启动
  - 其参数类型是：<<<dim3 gridDim, dim3 blockDim, size_t bytesSharedMemorySize, cudaStream_t stream>>>
    - dim3有构造函数dim3(int x, int y=1, int z=1)
      - 因此当直接赋值为int时，实则定义了dim.x = value, dim.y = 1, dim.z = 1
  - 其中gridDim, blockDim, bytesSharedMemory, stream是线程layout参数
    - 如果指定了stream，则把核函数加入到stream中异步执行
  - pdata和ndata则是核函数的函数调用参数
  - 函数调用参数必须传值，不能传引用等。参数可以是类类型等
- 核函数的执行无论stream是否为nullptr，都将是异步执行
  - 因此在核函数中进行printf操作，你必须进行等待，例如cudaDeviceSynchronize、或者cudaStreamSynchronize，否则你将无法看到打印的信息
__device__标记，设备调用的函数
- 调用方必须是device
__host__标记，主机调用函数
- 调用方必须是主机
也可以__device__ __host__两个标记同时有，表明该函数可以设备也可以主机
__constant__标记，定义常量内存
__shared__标记，定义共享内存
可以通过cudaPeekAtLastError/cudaGetLastError函数，捕获核函数是否出现错误或者异常

（五）核函数的线程总数

host调用核函数：function<<<gridDim, blockDim, sharedMemorySize, stream>>>(args…);

在核函数内部有四个非常有用的内置变量——threadIdx、blockIdx、blockDim和gridDim。线程layout主要用到blockDim、gridDim。

首先，我们可以先不严谨地认为GPU是一个立方体，这个立方体有很多小方块：

之后其中的每一个小块都是一个thread，为了方便讨论，我们只考虑2D：

我们关心的是某一个thread的位置，比如图上的黄色方块：

其位置为(blockIdx.x, blockIdx.y, threadIdx.x, threadIdx.y) =(1, 0, 1, 1)，如果把这个2D转换为1D，这格黄色的thread的1D位置为13。

我们可以通过blockIdx索引到线程块，通过threadIdx索引到某个块内的线程，通过blockDim得到一个块内线程总数，通过gridDim得到一个格内块总数：

那么我们可以得到在一维的条件下，计算线程全局id公式为：

线程全局id = blockIdex.x * blockDim.x + threadIdx.x

其表示的含义是要求thread的1D idx，先得知道在第几个block里，再知道在这个block里的第几个thread。

在一维的条件下，计算核函数内的线程总数为：

核函数的线程总数 = gridDim.x * blockDim.x

那我们就可以暗中维度高低排序来看待这个信息：

之后计算的话就这么计算：（左乘右加）

（六）代码

kernel.cu：

#include <stdio.h>
#include <cuda_runtime.h>

__global__ void test_print_kernel(const float* pdata, int ndata){

    int idx = threadIdx.x + blockIdx.x * blockDim.x;
    /*    dims                 indexs
        gridDim.z            blockIdx.z
        gridDim.y            blockIdx.y
        gridDim.x            blockIdx.x
        blockDim.z           threadIdx.z
        blockDim.y           threadIdx.y
        blockDim.x           threadIdx.x

        Pseudo code:
        position = 0
        for i in 6:
            position *= dims[i]
            position += indexs[i]
    */
    printf("Element[%d] = %f, threadIdx.x=%d, blockIdx.x=%d, blockDim.x=%d\n", idx, pdata[idx], threadIdx.x, blockIdx.x, blockDim.x);
}

void test_print(const float* pdata, int ndata){

    // <<<gridDim, blockDim, bytes_of_shared_memory, stream>>>
    test_print_kernel<<<1, ndata, 0, nullptr>>>(pdata, ndata);

    // 在核函数执行结束后，通过cudaPeekAtLastError获取得到的代码，来知道是否出现错误
    // cudaPeekAtLastError和cudaGetLastError都可以获取得到错误代码
    // cudaGetLastError是获取错误代码并清除掉，也就是再一次执行cudaGetLastError获取的会是success
    // 而cudaPeekAtLastError是获取当前错误，但是再一次执行 cudaPeekAtLastError 或者 cudaGetLastError 拿到的还是那个错
    // cuda的错误会传递，如果这里出错了，不移除。那么后续的任意api的返回值都会是这个错误，都会失败
    cudaError_t code = cudaPeekAtLastError();
    if(code != cudaSuccess){    
        const char* err_name    = cudaGetErrorName(code);    
        const char* err_message = cudaGetErrorString(code);  
        printf("kernel error %s:%d  test_print_kernel failed. \n  code = %s, message = %s\n", __FILE__, __LINE__, err_name, err_message);   
    }
}

main.cpp：


#include <cuda_runtime.h>
#include <stdio.h>

#define checkRuntime(op)  __check_cuda_runtime((op), #op, __FILE__, __LINE__)

bool __check_cuda_runtime(cudaError_t code, const char* op, const char* file, int line){
    if(code != cudaSuccess){    
        const char* err_name = cudaGetErrorName(code);    
        const char* err_message = cudaGetErrorString(code);  
        printf("runtime error %s:%d  %s failed. \n  code = %s, message = %s\n", file, line, op, err_name, err_message);   
        return false;
    }
    return true;
}

void test_print(const float* pdata, int ndata);

int main(){
    float* parray_host = nullptr;
    float* parray_device = nullptr;
    int narray = 10;
    int array_bytes = sizeof(float) * narray;

    parray_host = new float[narray];
    checkRuntime(cudaMalloc(&parray_device, array_bytes));

    for(int i = 0; i < narray; ++i)
        parray_host[i] = i;
    
    checkRuntime(cudaMemcpy(parray_device, parray_host, array_bytes, cudaMemcpyHostToDevice));
    test_print(parray_device, narray);
    checkRuntime(cudaDeviceSynchronize());

    checkRuntime(cudaFree(parray_device));
    delete[] parray_host;
    return 0;
}

0x08 共享内存

共享内存特性：

共享内存的主要特点在于“共享”，也即同一个线程块中的所有线程都可以对这一块存储进行读写操作，所以“共享”是针对同一个线程块所有线程而言的。一旦共享内存被定义并指定大小，系统将给所有线程块都分配相同大小的共享内存。

共享内存更靠近计算单元，所以访问速度更快。
共享内存通常可以作为访问全局内存的缓存使用。
可以利用共享内存实现线程间的通信。
通常与__syncthreads同时出现，这个函数是同步block内的所有线程，全部执行到这一行才往下走。
使用方式，通常是在线程id为0的时候从global memory取值，然后syncthreads，然后再使用。

共享内存在CUDA核函数中定义，通常有两种方式：静态方式、动态方式。

静态方式定义：这种定义的特点是定义的同时指定大小，并且定义内存大小通常地址是不一样的。

const size_t static_shared_memory_num_element = 6 * 1024; // 6KB
__shared__ char static_shared_memory[static_shared_memory_num_element]; 
__shared__ char static_shared_memory2[2];

动态方式定义：这种定义的特点是定义的时候不指定大小，在调用核函数的时候将共享内存大小以输入参数的形式传入。动态定义共享内存时，是使用核函数进行传入数值的，传入的数值必须以字节byte为单位。如果定义的共享内存不是byte类型，数值必须乘以类型所占用的字节数。

// 表示的是传入12个字节大小的共享内存
demo1_kernel<<<1, 1, 12, nullptr>>>();

如何知道本机共享内存的大小：

sharedMemPerBlock 指示了block中最大可用的共享内存，共享内存使得block内存内的threads可以相互通信。

共享内存是片上内存，更靠近计算单元，因此比globalMem速度更快，通常可以充当缓存使用。因此在运算的时候，可以数据先读入到sharedMem，做各类计算时，使用sharedMem而非globalMem。

动态共享内存的使用注意：
- dynamic_shared_memory变量必须使用extern shared开头
- 并且定义为不确定大小的数组[]
- 变量放在函数外面和里面都一样
- 其指针由cuda调度器执行时赋值
- 动态共享内存无论定义多少个，地址都一样
静态共享内存的使用注意：
- 不加extern，以shared开头
- 定义时需要明确数组的大小
- 静态分配的地址比动态分配的地址低
- 静态共享内存定义几个地址随之叠加
如果配置的各类共享内存总和大于sharedMemPerBlock，则核函数执行出错，Invalid argument。
- 不同类型的静态共享变量定义，其内存划分并不一定是连续的。
- 中间会有内存对齐策略，使得第一个和第二个变量之间可能存在空隙。
- 因此你的变量之间如果存在空隙，可能小于全部大小的共享内存就会报错。

代码：

shared-memory.cu：

#include <cuda_runtime.h>
#include <stdio.h>

//demo1 //
/* 
demo1 主要为了展示查看静态和动态共享变量的地址
 */
const size_t static_shared_memory_num_element = 6 * 1024; // 6KB
__shared__ char static_shared_memory[static_shared_memory_num_element]; 
__shared__ char static_shared_memory2[2]; 

__global__ void demo1_kernel(){
    extern __shared__ char dynamic_shared_memory[];      // 静态共享变量和动态共享变量在kernel函数内/外定义都行，没有限制
    extern __shared__ char dynamic_shared_memory2[];
    printf("static_shared_memory = %p\n",   static_shared_memory);   // 静态共享变量，定义几个地址随之叠加
    printf("static_shared_memory2 = %p\n",  static_shared_memory2); 
    printf("dynamic_shared_memory = %p\n",  dynamic_shared_memory);  // 动态共享变量，无论定义多少个，地址都一样
    printf("dynamic_shared_memory2 = %p\n", dynamic_shared_memory2); 

    if(blockIdx.x == 0 && threadIdx.x == 0) // 第一个thread
        printf("Run kernel.\n");
}

/demo2//
/* 
demo2 主要是为了演示的是如何给 共享变量进行赋值
 */
// 定义共享变量，但是不能给初始值，必须由线程或者其他方式赋值
__shared__ int shared_value1;

__global__ void demo2_kernel(){
    
    __shared__ int shared_value2;
    if(threadIdx.x == 0){

        // 在线程索引为0的时候，为shared value赋初始值
        if(blockIdx.x == 0){
            shared_value1 = 123;
            shared_value2 = 55;
        }else{
            shared_value1 = 331;
            shared_value2 = 8;
        }
    }

    // 等待block内的所有线程执行到这一步
    __syncthreads();
    
    printf("%d.%d. shared_value1 = %d[%p], shared_value2 = %d[%p]\n", 
        blockIdx.x, threadIdx.x,
        shared_value1, &shared_value1, 
        shared_value2, &shared_value2
    );
}

void launch(){
    
    demo1_kernel<<<1, 1, 12, nullptr>>>();
    demo2_kernel<<<2, 5, 0, nullptr>>>();
}

main.cpp：

#include <cuda_runtime.h>
#include <stdio.h>

#define checkRuntime(op)  __check_cuda_runtime((op), #op, __FILE__, __LINE__)

bool __check_cuda_runtime(cudaError_t code, const char* op, const char* file, int line){
    if(code != cudaSuccess){    
        const char* err_name = cudaGetErrorName(code);    
        const char* err_message = cudaGetErrorString(code);  
        printf("runtime error %s:%d  %s failed. \n  code = %s, message = %s\n", file, line, op, err_name, err_message);   
        return false;
    }
    return true;
}

void launch();

int main(){

    cudaDeviceProp prop;
    checkRuntime(cudaGetDeviceProperties(&prop, 0));
    printf("prop.sharedMemPerBlock = %.2f KB\n", prop.sharedMemPerBlock / 1024.0f);

    launch();
    checkRuntime(cudaPeekAtLastError());
    checkRuntime(cudaDeviceSynchronize());
    printf("done\n");
    return 0;
}

0x09 Warpaffine

主要解决图像的缩放喝平移来处理目标检测中常见的预处理行为：

warpaffine是对图像做平移缩放旋转进行综合统一描述的方法，同时也是一个很容易实现cuda并行加速的算法。
在深度学习领域通常需要做预处理，比如CopyMakeBorder，RGB->BGR，减去均值除以标准差，BGRBGRBGR -> BBBGGGRRR。
如果使用cuda进行并行加速实现，那么整个预处理都进行统一，并且性能也很好。
由于warpaffine是标准的矩阵映射坐标，并且可逆，所以逆变换就是其变换矩阵的逆矩阵。

warpaffine本质：对原图的坐标使用矩阵进行变换，之后生成新的图。

（一）点的缩放平移旋转

点的变换有：缩放、平移、旋转。假设我们现在要对点P(x,y)进行旋转θ度，缩放scale倍，平移ox、oy，我们可以这么去看：

先进行旋转变换，矩阵可以这么表示：

其详细的坐标变换：

又因为opencv的图像坐标，原点在左上角，y+向下，因此旋转变换在矩阵上的图像是：

之后再看看缩放变换：

那么缩放+旋转变换：

平移变换：

这种变换是加法运算，无法合并为一个矩阵，这个时候我们就需要齐次坐标定义：用(x, y, w)，来表示点P(x/w,y/w,1)

我们先将平移变换转换为这个矩阵：

这个时候缩放+旋转+平移，要注意顺序：

旋转矩阵的逆变换，可以说是其转置:

反变换，对变换矩阵求逆即可：

对于平常的目标检测推理而言，通常需要图像等比缩放并且居中，那么可以分为以下三步：

因此根据上面的公式，我们可以看到变换可以表示为：P'=TOSP。

由于三个变换比较简单，我们可以直接写出矩阵：

那么逆变换：

（二）线性插值

线性插值是指使用连接两个已知量的直线来确定在这两个已知量之间的一个未知量值的方法。

假设我们已经知道坐标（x0，y0）以及（x1，y1），我们要得到[x0,x1]区间内某一位置x再直线上的值，如图：

那么我们就可以通过公式求出这条两点式方程：

假设方程两边的值为a，那么这个值就是差值系数从x0到x1距离的比值，由于x值抑制，所以我们可以再公式中得到：

同理：

那么这样在代数上就可以表示为：

或者是

这样我们就可以通过a来得到未知值y。

（三）双线性差值

假设我们要求这个紫色点的像素值：

要求这个紫色点的颜色我们就需要知道它周围四个点颜色的加权和，而每个点的权重，则是其对面矩形区域的面积，占总面积的比例。

其算法原理为，我们先得到我们周围的四个点，使用保留最近的最大整数值并且小于我们当前的坐标x，求得x1，x2则使用+1的形式进行求解，那么y1、y2也是同理。这样我们就可以得出这个周围框框的范围，之后我们再使用四个点的像素值，去乘以它们对角的面积，也就是上面那张图那样的求法，就可以得出相应的值。

那么使用python实现的话可以这么写：

def pyWarpAffine(image, M, dst_size, constant=(0, 0, 0)):
    
    # 注意输入的M矩阵格式，是Origin->Dst
    # 而这里需要的是Dst->Origin，所以要取逆矩阵
    M = cv2.invertAffineTransform(M)
    constant = np.array(constant)
    ih, iw   = image.shape[:2]
    dw, dh   = dst_size
    dst      = np.full((dh, dw, 3), constant, dtype=np.uint8)
    irange   = lambda p: p[0] >= 0 and p[0] < iw and p[1] >= 0 and p[1] < ih
    
    for y in range(dh):
        for x in range(dw):
            
            homogeneous = np.array([[x, y, 1]]).T
            ox, oy = M @ homogeneous
            
            low_ox = int(np.floor(ox))
            low_oy = int(np.floor(oy))
            high_ox = low_ox + 1
            high_oy = low_oy + 1
            
            # p0     p1
            #      o
            # p2     p3
            
            pos = ox - low_ox, oy - low_oy
            # 求出各个区域的面积
            p0_area = (1 - pos[0]) * (1 - pos[1])
            p1_area = pos[0] * (1 - pos[1])
            p2_area = (1 - pos[0]) * pos[1]
            p3_area = pos[0] * pos[1]
            
            p0 = low_ox, low_oy
            p1 = high_ox, low_oy
            p2 = low_ox, high_oy
            p3 = high_ox, high_oy
            p0_value = image[p0[1], p0[0]] if irange(p0) else constant
            p1_value = image[p1[1], p1[0]] if irange(p1) else constant
            p2_value = image[p2[1], p2[0]] if irange(p2) else constant
            p3_value = image[p3[1], p3[0]] if irange(p3) else constant
            dst[y, x] = p0_area * p0_value + p1_area * p1_value + p2_area * p2_value + p3_area * p3_value
    return dst

看看cpp版本：

__global__ void warp_affine_bilinear_kernel(
    uint8_t* src, int src_line_size, int src_width, int src_height, 
    uint8_t* dst, int dst_line_size, int dst_width, int dst_height, 
	uint8_t fill_value, AffineMatrix matrix
){   
    int dx = blockDim.x * blockIdx.x + threadIdx.x; 
    int dy = blockDim.y * blockIdx.y + threadIdx.y;
    if (dx >= dst_width || dy >= dst_height)  return;

    float c0 = fill_value, c1 = fill_value, c2 = fill_value;
    float src_x = 0; float src_y = 0;
    affine_project(matrix.d2i, dx, dy, &src_x, &src_y);

    if(src_x < -1 || src_x >= src_width || src_y < -1 || src_y >= src_height){
        // out of range
        // src_x < -1时，其高位high_x < 0，超出范围
        // src_x >= -1时，其高位high_x >= 0，存在取值
    }else{
		// 临近最大的整数值，但不大于本身
        int y_low = floorf(src_y);
        int x_low = floorf(src_x);
        int y_high = y_low + 1;
        int x_high = x_low + 1;
		// 这个是要扩充的像素值 代表RGB
        uint8_t const_values[] = {fill_value, fill_value, fill_value};
        float ly    = src_y - y_low;
        float lx    = src_x - x_low;
        float hy    = 1 - ly;
        float hx    = 1 - lx;
        float w1    = hy * hx, w2 = hy * lx, w3 = ly * hx, w4 = ly * lx;
		// 这个constvalue其实就是对图像扩充时所需要用到的像素值
        uint8_t* v1 = const_values;
        uint8_t* v2 = const_values;
        uint8_t* v3 = const_values;
        uint8_t* v4 = const_values;
        if(y_low >= 0){
			//这里要确定当y满足在图像框中时，x的两个坐标坐标是否超出图像范围
            if (x_low >= 0)
				// 获取该像素的地址，这个时候使用v[1]、v[2]、v[3]就可以取到RGB
                v1 = src + y_low * src_line_size + x_low * 3;

            if (x_high < src_width)
                v2 = src + y_low * src_line_size + x_high * 3;
        }
        
        if(y_high < src_height){
            if (x_low >= 0)
                v3 = src + y_high * src_line_size + x_low * 3;

            if (x_high < src_width)
                v4 = src + y_high * src_line_size + x_high * 3;
        }
        
		// 对临近的四个点进行加权求和 得到新的rgb值
        c0 = floorf(w1 * v1[0] + w2 * v2[0] + w3 * v3[0] + w4 * v4[0] + 0.5f);
        c1 = floorf(w1 * v1[1] + w2 * v2[1] + w3 * v3[1] + w4 * v4[1] + 0.5f);
        c2 = floorf(w1 * v1[2] + w2 * v2[2] + w3 * v3[2] + w4 * v4[2] + 0.5f);
    }

    uint8_t* pdst = dst + dy * dst_line_size + dx * 3;
    pdst[0] = c0; pdst[1] = c1; pdst[2] = c2;
}

（五）代码

main.cpp：


#include <cuda_runtime.h>
#include <opencv2/opencv.hpp>
#include <stdio.h>

using namespace cv;

#define min(a, b)  ((a) < (b) ? (a) : (b))
#define checkRuntime(op)  __check_cuda_runtime((op), #op, __FILE__, __LINE__)

bool __check_cuda_runtime(cudaError_t code, const char* op, const char* file, int line){
    if(code != cudaSuccess){
        const char* err_name = cudaGetErrorName(code);    
        const char* err_message = cudaGetErrorString(code);  
        printf("runtime error %s:%d  %s failed. \n  code = %s, message = %s\n", file, line, op, err_name, err_message);   
        return false;
    }
    return true;
}

void warp_affine_bilinear( // 声明
    uint8_t* src, int src_line_size, int src_width, int src_height, 
    uint8_t* dst, int dst_line_size, int dst_width, int dst_height, 
	uint8_t fill_value
);

Mat warpaffine_to_center_align(const Mat& image, const Size& size){  
    
    Mat output(size, CV_8UC3);
    // 声明两个指针
    uint8_t* psrc_device = nullptr;
    uint8_t* pdst_device = nullptr;
    size_t src_size = image.cols * image.rows * 3;
    size_t dst_size = size.width * size.height * 3;
	
    // 在GPU上声明两个地址，第一个存放原始图像，第二个存放处理后的图像，之后将生成的图像再返回给CPU
    checkRuntime(cudaMalloc(&psrc_device, src_size)); // 在GPU上开辟两块空间
    checkRuntime(cudaMalloc(&pdst_device, dst_size));
    // 搬运数据到GPU上
    checkRuntime(cudaMemcpy(psrc_device, image.data, src_size, cudaMemcpyHostToDevice)); 
    
    warp_affine_bilinear(
        psrc_device, image.cols * 3, image.cols, image.rows,
        pdst_device, size.width * 3, size.width, size.height,
        114
    );

    // 检查核函数执行是否存在错误
    checkRuntime(cudaPeekAtLastError());
    // 将预处理完的数据搬运回来
    checkRuntime(cudaMemcpy(output.data, pdst_device, dst_size, cudaMemcpyDeviceToHost)); 	  // 要记得释放GPU上的内存
    checkRuntime(cudaFree(psrc_device));
    checkRuntime(cudaFree(pdst_device));
    return output;
}

int main(){ 
    Mat image = imread("xxx.jpg");
    Mat output = warpaffine_to_center_align(image, Size(640, 640));
    imwrite("output.jpg", output);
    printf("Done. save to output.jpg\n");
    return 0;
}

再看看处理函数affine.cu：

#include <cuda_runtime.h>

#define min(a, b)  ((a) < (b) ? (a) : (b))
#define num_threads   512

typedef unsigned char uint8_t;

struct Size{
    int width = 0, height = 0;

	//恢复默认构造
    Size() = default;
    Size(int w, int h)
    :width(w), height(h){}
};

// 计算仿射变换矩阵
// 计算的矩阵是居中缩放
struct AffineMatrix{
    
    float i2d[6];       // image to dst(network), 2x3 matrix
    float d2i[6];       // dst to image, 2x3 matrix

    // 这里其实是求解imat的逆矩阵，由于这个3x3矩阵的第三行是确定的0, 0, 1，因此可以简写如下
    void invertAffineTransform(float imat[6], float omat[6]){
        float i00 = imat[0];  float i01 = imat[1];  float i02 = imat[2];
        float i10 = imat[3];  float i11 = imat[4];  float i12 = imat[5];

        // 计算行列式
        float D = i00 * i11 - i01 * i10;
        D = D != 0 ? 1.0 / D : 0;

        // 计算剩余的伴随矩阵除以行列式
        float A11 = i11 * D;
        float A22 = i00 * D;
        float A12 = -i01 * D;
        float A21 = -i10 * D;
        float b1 = -A11 * i02 - A12 * i12;
        float b2 = -A21 * i02 - A22 * i12;
        omat[0] = A11;  omat[1] = A12;  omat[2] = b1;
        omat[3] = A21;  omat[4] = A22;  omat[5] = b2;
    }

    void compute(const Size& from, const Size& to){
        float scale_x = to.width / (float)from.width;
        float scale_y = to.height / (float)from.height;

        // 这里取min的理由是
        // 1. M矩阵是 from * M = to的方式进行映射，因此scale的分母一定是from
        // 2. 取最小，即根据宽高比，算出最小的比例，如果取最大，则势必有一部分超出图像范围而被裁剪掉，这不是我们要的
        float scale = min(scale_x, scale_y); 
        /**
        这里的仿射变换矩阵实质上是2x3的矩阵，具体实现是
        scale, 0, -scale * from.width * 0.5 + to.width * 0.5
        0, scale, -scale * from.height * 0.5 + to.height * 0.5
        
        这里可以想象成，是经历过缩放、平移、平移三次变换后的组合，M = TPS
        例如第一个S矩阵，定义为把输入的from图像，等比缩放scale倍，到to尺度下
        S = [
        scale,     0,      0
        0,     scale,      0
        0,         0,      1
        ]
        
        P矩阵定义为第一次平移变换矩阵，将图像的原点，从左上角，移动到缩放(scale)后图像的中心上
        P = [
        1,        0,      -scale * from.width * 0.5
        0,        1,      -scale * from.height * 0.5
        0,        0,                1
        ]

        T矩阵定义为第二次平移变换矩阵，将图像从原点移动到目标（to）图的中心上
        T = [
        1,        0,      to.width * 0.5,
        0,        1,      to.height * 0.5,
        0,        0,            1
        ]

        通过将3个矩阵顺序乘起来，即可得到下面的表达式：
        M = [
        scale,    0,     -scale * from.width * 0.5 + to.width * 0.5
        0,     scale,    -scale * from.height * 0.5 + to.height * 0.5
        0,        0,                     1
        ]
        去掉第三行就得到opencv需要的输入2x3矩阵
        **/

        i2d[0] = scale;  i2d[1] = 0;  i2d[2] = 
            -scale * from.width  * 0.5  + to.width * 0.5 + scale * 0.5 - 0.5;

        i2d[3] = 0;  i2d[4] = scale;  i2d[5] = 
            -scale * from.height * 0.5 + to.height * 0.5 + scale * 0.5 - 0.5;

        invertAffineTransform(i2d, d2i);
    }
};

__device__ void affine_project(float* matrix, int x, int y, float* proj_x, float* proj_y){

    // matrix
    // m0, m1, m2
    // m3, m4, m5
    *proj_x = matrix[0] * x + matrix[1] * y + matrix[2];
    *proj_y = matrix[3] * x + matrix[4] * y + matrix[5];
}

__global__ void warp_affine_bilinear_kernel(
    uint8_t* src, int src_line_size, int src_width, int src_height, 
    uint8_t* dst, int dst_line_size, int dst_width, int dst_height, 
	uint8_t fill_value, AffineMatrix matrix
){
    int dx = blockDim.x * blockIdx.x + threadIdx.x; 
    int dy = blockDim.y * blockIdx.y + threadIdx.y;
    if (dx >= dst_width || dy >= dst_height)  return;

    float c0 = fill_value, c1 = fill_value, c2 = fill_value;
    float src_x = 0; float src_y = 0;
    affine_project(matrix.d2i, dx, dy, &src_x, &src_y);

    if(src_x < -1 || src_x >= src_width || src_y < -1 || src_y >= src_height){
        // out of range
        // src_x < -1时，其高位high_x < 0，超出范围
        // src_x >= -1时，其高位high_x >= 0，存在取值
    }else{
		// 临近最大的整数值，但不大于本身
        int y_low = floorf(src_y);
        int x_low = floorf(src_x);
        int y_high = y_low + 1;
        int x_high = x_low + 1;

        uint8_t const_values[] = {fill_value, fill_value, fill_value};
        float ly    = src_y - y_low;
        float lx    = src_x - x_low;
        float hy    = 1 - ly;
        float hx    = 1 - lx;
        float w1    = hy * hx, w2 = hy * lx, w3 = ly * hx, w4 = ly * lx;
		// 这个constvalue其实就是对图像扩充时所需要用到的像素值
        uint8_t* v1 = const_values;
        uint8_t* v2 = const_values;
        uint8_t* v3 = const_values;
        uint8_t* v4 = const_values;
        if(y_low >= 0){
			//这里要确定当y满足在图像框中时，x的两个坐标坐标是否超出图像范围
            if (x_low >= 0)
				// 获取该像素的地址，这个时候使用v[1]、v[2]、v[3]就可以取到RGB
                v1 = src + y_low * src_line_size + x_low * 3;

            if (x_high < src_width)
                v2 = src + y_low * src_line_size + x_high * 3;
        }
        
        if(y_high < src_height){
            if (x_low >= 0)
                v3 = src + y_high * src_line_size + x_low * 3;

            if (x_high < src_width)
                v4 = src + y_high * src_line_size + x_high * 3;
        }
        
		// 对临近的四个点进行加权求和 得到新的rgb值
        c0 = floorf(w1 * v1[0] + w2 * v2[0] + w3 * v3[0] + w4 * v4[0] + 0.5f);
        c1 = floorf(w1 * v1[1] + w2 * v2[1] + w3 * v3[1] + w4 * v4[1] + 0.5f);
        c2 = floorf(w1 * v1[2] + w2 * v2[2] + w3 * v3[2] + w4 * v4[2] + 0.5f);
    }

    uint8_t* pdst = dst + dy * dst_line_size + dx * 3;
    pdst[0] = c0; pdst[1] = c1; pdst[2] = c2;
}

void warp_affine_bilinear(
    uint8_t* src, int src_line_size, int src_width, int src_height, 
    uint8_t* dst, int dst_line_size, int dst_width, int dst_height, 
	uint8_t fill_value
){
    dim3 block_size(32, 32); // blocksize最大就是1024，这里用2d来看更好理解
    dim3 grid_size((dst_width + 31) / 32, (dst_height + 31) / 32);
    AffineMatrix affine;
    affine.compute(Size(src_width, src_height), Size(dst_width, dst_height));

    warp_affine_bilinear_kernel<<<grid_size, block_size, 0, nullptr>>>(
        src, src_line_size, src_width, src_height,
        dst, dst_line_size, dst_width, dst_height,
        fill_value, affine
    );
}

0x0A Yolov5加速后处理

对于模型推理后的后处理，可以直接使用cuda核函数进行解码，效率比较高。

那么加速处理的重点在哪：

CPU解码重点：

避免多余的计算，需要知道有些数学运算需要的时间远超过很多if，减少他们的计算次数就是性能的关键。

nms的实现是可以优化的，例如remove flag并且预先分配内存，reserve对输出分配内存。

GPU解码重点：

表示输出数量不确定的数组，用[count, box1, box2, box3]的方式，此时需要有最大数量限制。

通过atomicAdd实现数组元素的加入，并返回索引。

一样的，不必要的计算，尽量省掉。

在main.cpp中，我们可以通过如下来减少计算的量：

vector<Box> cpu_decode(float* predict, int rows, int cols, float confidence_threshold = 0.25f, float nms_threshold = 0.45f){
    
    vector<Box> boxes;
    int num_classes = cols - 5;
    for(int i = 0; i < rows; ++i){
        float* pitem = predict + i * cols;
        float objness = pitem[4];
		// 这么处理可以节省大量时间
        if(objness < confidence_threshold)
            continue;

        float* pclass = pitem + 5;
        int label     = std::max_element(pclass, pclass + num_classes) - pclass;
        float prob    = pclass[label];
        float confidence = prob * objness;
		// 这么处理可以节省大量时间
        if(confidence < confidence_threshold)
            continue;

		// 滤掉不符合的数据 再进行储存
        float cx     = pitem[0];
        float cy     = pitem[1];
        float width  = pitem[2];
        float height = pitem[3];
        float left   = cx - width * 0.5;
        float top    = cy - height * 0.5;
        float right  = cx + width * 0.5;
        float bottom = cy + height * 0.5;
        boxes.emplace_back(left, top, right, bottom, confidence, (float)label);
    }

    std::sort(boxes.begin(), boxes.end(), [](Box& a, Box& b){return a.confidence > b.confidence;});
	// 使用vector<bool>来标记是否需要删除box
    std::vector<bool> remove_flags(boxes.size());
    std::vector<Box> box_result;
	// 预分配内存
    box_result.reserve(boxes.size());

	// iou计算函数 计算两个box之间的iou
    auto iou = [](const Box& a, const Box& b){
        float cross_left   = std::max(a.left, b.left);
        float cross_top    = std::max(a.top, b.top);
        float cross_right  = std::min(a.right, b.right);
        float cross_bottom = std::min(a.bottom, b.bottom);
		// 计算重叠面积
        float cross_area = std::max(0.0f, cross_right - cross_left) * std::max(0.0f, cross_bottom - cross_top);
		// A面积+B面积
        float union_area = std::max(0.0f, a.right - a.left) * std::max(0.0f, a.bottom - a.top) 
                         + std::max(0.0f, b.right - b.left) * std::max(0.0f, b.bottom - b.top) - cross_area;
        if(cross_area == 0 || union_area == 0) return 0.0f;
        return cross_area / union_area;
    };

    for(int i = 0; i < boxes.size(); ++i){
		// 判断是否需要执行iou大小的比对
        if(remove_flags[i]) continue;
		// 第一个box肯定不会remove
        auto& ibox = boxes[i];
		// 插入
        box_result.emplace_back(ibox);
        for(int j = i + 1; j < boxes.size(); ++j){
			// 判断是否为移除box
            if(remove_flags[j]) continue;

            auto& jbox = boxes[j];
			// 如果它们的标签相同，那么比较iou，移除其中一个
            if(ibox.label == jbox.label){
                // class matched
                if(iou(ibox, jbox) >= nms_threshold)
                    remove_flags[j] = true;
            }
        }
    }
    return box_result;
}

在GPU中的优化：

void decode_kernel_invoker(
    float* predict, int num_bboxes, int num_classes, float confidence_threshold, 
    float nms_threshold, float* invert_affine_matrix, float* parray, int max_objects, int NUM_BOX_ELEMENT, cudaStream_t stream);

vector<Box> gpu_decode(float* predict, int rows, int cols, float confidence_threshold = 0.25f, float nms_threshold = 0.45f){
    
    vector<Box> box_result;
    cudaStream_t stream = nullptr;
    checkRuntime(cudaStreamCreate(&stream));

    float* predict_device = nullptr;
    float* output_device = nullptr;
    float* output_host = nullptr;
	// 先申请内存 实际不可能有这么大
    int max_objects = 1000;
    int NUM_BOX_ELEMENT = 7;  // left, top, right, bottom, confidence, class, keepflag
    checkRuntime(cudaMalloc(&predict_device, rows * cols * sizeof(float)));
    checkRuntime(cudaMalloc(&output_device, sizeof(float) + max_objects * NUM_BOX_ELEMENT * sizeof(float)));
    checkRuntime(cudaMallocHost(&output_host, sizeof(float) + max_objects * NUM_BOX_ELEMENT * sizeof(float)));

    checkRuntime(cudaMemcpyAsync(predict_device, predict, rows * cols * sizeof(float), cudaMemcpyHostToDevice, stream));
    decode_kernel_invoker(
        predict_device, rows, cols - 5, confidence_threshold, 
        nms_threshold, nullptr, output_device, max_objects, NUM_BOX_ELEMENT, stream
    );
    checkRuntime(cudaMemcpyAsync(output_host, output_device, 
        sizeof(int) + max_objects * NUM_BOX_ELEMENT * sizeof(float), 
        cudaMemcpyDeviceToHost, stream
    ));
    checkRuntime(cudaStreamSynchronize(stream));

    int num_boxes = min((int)output_host[0], max_objects);
    for(int i = 0; i < num_boxes; ++i){
        float* ptr = output_host + 1 + NUM_BOX_ELEMENT * i;
        int keep_flag = ptr[6];
        if(keep_flag){
            box_result.emplace_back(
                ptr[0], ptr[1], ptr[2], ptr[3], ptr[4], (int)ptr[5]
            );
        }
    }
	// 要记得释放内存
    checkRuntime(cudaStreamDestroy(stream));
    checkRuntime(cudaFree(predict_device));
    checkRuntime(cudaFree(output_device));
    checkRuntime(cudaFreeHost(output_host));
    return box_result;
}

0x0B Thrust and Error

在这里只是补充了几个小知识点：

thrust是CUDA开发的，可以理解为cpp中的stl库，因为通常没用到thrust，所以就在这里当一个小知识点：

对于thrust中的lambda表达式，需要增加device标记表明函数可以被核函数调用，此时需要在makefile中增加--extended-lambda标记
由于使用到了device vector，因此编译环境需要修改为nvcc编译，因此main.cpp改成了main.cu
内存的复制和分配，被cuda封装

其使用如下：

#include <stdio.h>
#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <thrust/sort.h>
#include <iostream>
using namespace std;

__host__ __device__
int sort_func(int a, int b){
    return a > b;
}

int main(){

    int data[] = {5, 3, 1, 5, 2, 0};
    int ndata  = sizeof(data) / sizeof(data[0]);
    thrust::host_vector<int> array1(data, data + ndata);
    thrust::sort(array1.begin(), array1.end(), sort_func);

    thrust::device_vector<int> array2 = thrust::host_vector<int>(data, data + ndata);
    thrust::sort(array2.begin(), array2.end(), []__device__(int a, int b){return a < b;});

    printf("array1------------------------\n");
    for(int i = 0; i < array1.size(); ++i)
        cout << array1[i] << endl;

    printf("array2------------------------\n");
    for(int i = 0; i < array2.size(); ++i)
        cout << array2[i] << endl;
    return 0;
}

那么关于cuda中的错误处理：

若cuda核函数出错，由于他是异步的，立即执行cudaPeekAtLastError只会拿到对输入参数校验是否正确的状态，而不会拿到核函数是否执行正确的状态。因此需要等待核函数执行完毕后，才真的知道当前核函数是否出错，一般通过设备同步或者流同步进行等待。

错误可以分为可恢复和不可恢复两种：

可恢复：

参数配置错误等，例如block越界（一般最大值是1024），shared memory大小超出范围（一般是48KB），也就是说当你调用核函数时，线程的设置设置错了，这是可以恢复的，会在下一次调用的时候被恢复。
通过cudaGetLastError可以获取错误代码，同时把当前状态恢复为success。
该错误在调用核函数后可以立即通过cudaGetLastError/cudaPeekAtLastError拿到。
该错误在下一个函数调用的时候会覆盖。

不可恢复：

核函数执行错误，例如访问越界等等异常
该错误则会传递到之后的所有cuda操作上
错误状态通常需要等到核函数执行完毕才能够拿到，也就是有可能在后续的任何流程中突然异常（因为是异步的）

例如：

#include <cuda_runtime.h>
#include <stdio.h>
#include <iostream>
using namespace std;

__global__ void func(float* ptr){

    int pos = blockIdx.x * blockDim.x + threadIdx.x;
    if(pos == 999){
        ptr[999] = 5;
    }
}

int main(){

    float* ptr = nullptr;

    // 因为核函数是异步的，因此不会立即检查到他是否存在异常
    func<<<100, 10>>>(ptr);
    //func<<<100, 1050>>>(ptr);  //这个错误是可以恢复的
    auto code1 = cudaPeekAtLastError();
    cout << cudaGetErrorString(code1) << endl;

    // 对当前设备的核函数进行同步，等待执行完毕，可以发现过程是否存在异常
    auto code2 = cudaDeviceSynchronize();
    cout << cudaGetErrorString(code2) << endl;

    // 异常会一直存在，以至于后续的函数都会失败
    float* new_ptr = nullptr;
    auto code3 = cudaMalloc(&new_ptr, 100);
    cout << cudaGetErrorString(code3) << endl;
    return 0;
}