# cudaNoKernel **Repository Path**: iam002/cuda-no-kernel ## Basic Information - **Project Name**: cudaNoKernel - **Description**: 本项目提供一系列模板函数，无需要手动写核函数而实现cuda编程 - **Primary Language**: Unknown - **License**: GPL-3.0 - **Default Branch**: master - **Homepage**: None - **GVP Project**: No ## Statistics - **Stars**: 0 - **Forks**: 0 - **Created**: 2024-10-14 - **Last Updated**: 2024-10-20 ## Categories & Tags **Categories**: Uncategorized **Tags**: cuda ## README # CUDA NO KERNEL ## 介绍 CUDA 核函数的代码编写往往存在大量的重复冗余工作，如： 1. 每次调用需要手动设置的 block size 和 grid size 2. 核函数内部需要手动检查和处理线程索引越界的情况 3. 不同数据类型（float 或 double）实现相同的运算或功能 4. 实现一些常见的计算任务（如逐元素运算，规约操作、扫描操作等） 5. $\cdots$ 为此，本项目通过**模板编程**来编写通用的核函数，并将模板核函数进一步封装为普通的 C++ 函数。使用者只需要实现核心的并行计算功能模块，无需过度关注 CUDA 核函数的编写细节，即可实现 CUDA 并行功能。 > Tips: 我的实现不是最优版本，仅供参考。 ## 文件说明 - [cmake](cmake) 存放cmake宏函数定义文件，包含第三方库路径配置和编译辅助宏或函数定义 - [common](common) 存放源代码（非cuda） - [cudalib](cudalib) 存放 cuda 源代码 - [doc](doc) 存放文档 - [test](test) 存放测试代码 - [CMakeLists.txt](CMakeLists.txt) 项目主CMake文件 - [README.md](README.md) 项目说明文档 ## 编译与使用 ### 编译前的准备 1. 安装编译器：目前只在 Linux gcc 10.5 编译过，需要 C++11 标准 2. 安装 CMake（版本不低于 3.12） 3. 配置好你的 CUDA 环境 4. 推荐安装 VSCode，并使用 C/C++ Extension Pack 插件包进行开发 5. 根据需要修改 [cmake/thirdlibs/cuda_config.cmake](cmake/thirdlibs/cuda_config.cmake) 中的 CUDA 本地安装路径 `CUDA_PREFIX` ### 编译步骤 ![compile](img/compile.gif) ## 并行应用场景 ### 逐元素运算 elementwise #### 应用场景 最常见，逐元素求两个或两个以上的输入数组的和、差，乘积等等，然后返回一个与输入一样长的结果数组，当然我们的返回也可以是任意数量的输出数组，前者我称为 elementwise MISO（逐元素多入单出），后者称为 elementwise MIMO（逐元素多入多出）。 > 参考： > - https://zhuanlan.zhihu.com/p/591058808 > - https://github.com/BBuf/how-to-optim-algorithm-in-cuda/tree/master/elementwise #### 特点 1. 数组长度：输入数组与输出数组长度一致，输入数组与输出数组之间的数据类型没有限制 2. IO配对关系：输入数组中某一个元素与输出数组中的具有相同索引的元素进行配对，即第`i`个输入元素的计算结果将保存在输出数组的第`i`个元素； 3. 在某些情况下，还需要数组的索引作为输入参数。 #### 使用实例 模板函数的定义在头文件 [cudalib/public/hs_cuda_elemwise.cuh](cudalib/public/hs_cuda_elemwise.cuh)，本实例的代码在 [test/cudaElemwise](test/cudaElemwise/) 和 [test/cudaElemwiseMIMO](test/cudaElemwiseMIMO/)。 > 本项目定义 `RasterData` 来负责管理内存的分配和释放（有CPU和CUDA两个版本），具体可参见文档 [doc/hs_rasterdata.md](doc/hs_rasterdata.md) 以及 [doc/hs_rasterview.md](doc/hs_rasterview.md)。 > > 在明确我们的计算任务需求后，我们需要准备的是： > - 定义最小的计算单元（功能模块/算子） > - 选择相应的应用场景中的合适模板函数 ##### 实例(1) 逐元素相加 1. 在 [test/cudaElemwise/fun.cu](test/cudaElemwise/fun.cu) 中 `#include "hs_cuda_elemwise.cuh"` 2. 定义并行计算中的最小计算单元(功能模块)： ```cpp #include "hs_cuda_elemwise.cuh" template struct AddFunctor { /* 最小计算单元: 输入 x、y, 返回 x + y */ __device__ T operator()(T x, T y) const { return x + y; } }; ``` > 说明：我们需要定义一个结构体并重载`()`函数，该函数的形参数等于最小计算单元的输入参数的个数，例如这里的两数之和，形参数为 2，同时该函数需要返回两数之和。 > > 注意：函数定义前缀必须包含 `__device__`，且必须进行 `const` 进行修饰。 3. 调用 `hs_cuda_elemwise.cuh` 中对应的模板函数，这里的输入有两个，输入有一个，选择 `hs::cuda::elemwise::Binary` ```cpp void elemAdd(const hs::raster::RasterData& host_mat_a , const hs::raster::RasterData& host_mat_b , hs::raster::RasterData& host_mat_res) { /* 定义GPU端上的栅格数据 */ hscuda::raster::RasterData device_mat_a; hscuda::raster::RasterData device_mat_b; hscuda::raster::RasterData device_mat_c; /* 将输入数据拷贝到GPU */ hscuda::raster::host2device(host_mat_a.view(), device_mat_a); hscuda::raster::host2device(host_mat_b.view(), device_mat_b); /* 输出数据分配显存 */ int size_x = device_mat_a.view().sizeX(); int size_y = device_mat_a.view().sizeY(); device_mat_c.init(size_x, size_y); /* 使用二入一出模板函数 */ int N = device_mat_a.view().numel(); hs::cuda::elemwise::Binary(N , device_mat_a.view().pData() , device_mat_b.view().pData() , AddFunctor() , device_mat_c.view().pData() ); cudaDeviceSynchronize(); /* GPU输出数据拷贝到主机端 */ device_mat_c.toHost(host_mat_res); } ``` > 说明：像 `hs::cuda::elemwise::Binary` 的模板函数的形参组成都是一样的：先是带处理数据的长度(这里的`N`)，然后是任意数量的输入数据指针(紧跟后面的`*.pdata()`)，接着是功能算子(这里的`AddFunctor`)，最后是任意数量的输出数据指针(最后面的`*.pdata()`)。其中输入和输出数据指针的顺序与功能算账重载`()`函数的形参**顺序一致**。 > > 除了这个模板函数，还有 `hs::cuda::elemwise::Unary` 单入单出，`hs::cuda::elemwise::Ternary` 三入单出，`hs::cuda::elemwise::Quaternary` 四入单出，更多的话可能就不支持，或者自己定义一个，直接看对应的头文件照猫画虎即可。 ##### 实例(2) 实现 meshgrid 本实例实现了类似 matlab 的 `meshgrid` 函数功能，例如 `meshgrid(x:5, y:4)` 的输出为： ``` mat xx: 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 mat yy: 0 0 0 0 0 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 ``` 与实例(1)不同的是，最小的计算单元需要依赖数组的索引作为参数，因此可以选用 `hs::cuda::elemwiseIndex::Unary` 里的模板函数，这里实际上只需要一个索引作为输入，但是模板函数中的输入数组的最小个数为1，这里把输出数组作为输入数组，但实际上不起作用哈。 1. 在 [test/cudaElemwise/fun.cu](test/cudaElemwise/fun.cu) 定义并行计算中的最小计算单元(功能模块)： ```cpp struct FunctorMeshgridX { public: FunctorMeshgridX(int size_x, int size_y) : m_size_x(size_x), m_size_y(size_y) { } /* 本实例中不需要使用到 tmp */ __device__ int operator()(int id, int tmp) const { int yid = id / m_size_x; int xid = id - yid * m_size_x; return xid; } private: int m_size_x; int m_size_y; }; struct FunctorMeshgridY { public: FunctorMeshgridY(int size_x, int size_y) : m_size_x(size_x), m_size_y(size_y) { } /* 本实例中不需要使用到 tmp */ __device__ int operator()(int id, int tmp) const { int yid = id / m_size_x; return yid; } private: int m_size_x; int m_size_y; }; ``` > 说明：这里规定，如果并行运算单元需要索引信息，默认功能算子的重载`()`函数的**第一个形参**必须是 `int` 来表示该元素对应的索引值，其余部分与 elemwise 的功能算子一致。 2. 进一步的封装 `hs::cuda::elemwiseIndex::Unary` ```cpp void meshgrid(int size_x, int size_y, hs::raster::RasterData& host_mat_xx, hs::raster::RasterData& host_mat_yy) { /* GPU上分配输出显存 */ hscuda::raster::RasterData device_mat_xx(size_x, size_y); hscuda::raster::RasterData device_mat_yy(size_x, size_y); int len = size_x * size_y; hs::cuda::elemwiseIndex::Unary(len , device_mat_xx.view().pData() , FunctorMeshgridX(size_x, size_y) , device_mat_xx.view().pData() ); cudaDeviceSynchronize(); hs::cuda::elemwiseIndex::Unary(len , device_mat_yy.view().pData() , FunctorMeshgridY(size_x, size_y) , device_mat_yy.view().pData() ); cudaDeviceSynchronize(); device_mat_xx.toHost(host_mat_xx); device_mat_yy.toHost(host_mat_yy); } ``` ##### 实例(3) MIMO 模式 MIMO 模式下的逐元素运算的命名空间为 `hs::cuda::elemwise::MIMO`，下面实现这样一个分解功能: $$ \begin{aligned} y_1 &= (x_1 + x_2)/2 \\ y_2 &= (x_1 - x_2)/2 \\ \end{aligned} $$ 1. 在 [test/cudaElemwiseMIMO/fun.cu](test/cudaElemwiseMIMO/fun.cu) 定义并行计算中的最小计算单元(功能模块)： ```cpp #include "hs_cuda_elemwise.cuh" template struct OddEvenFunctor { __device__ void operator()(T x, T y, T & z1, T & z2) const { z1 = (x + y) / 2.0; z2 = (x - y) / 2.0; } }; ``` > 说明：之前的例子，只有一个输出值，可以使用 `return` 返回。如果存在多个输入，规定功能算子的重载`()`函数不再有返回值，输出以形参的形式表现在函数定义上，且所有输入形参排在所有输出形参的前面，形参的顺序与模板函数的数据指针顺序一致。 2. 这是一个二输入、二输出的计算单元，使用模板函数 `hs::cuda::elemwise::mimo::Op2in2out` 并做进一步的分装 ```cpp void testOddEven(const hs::raster::RasterData& host_mat_a , const hs::raster::RasterData& host_mat_b , hs::raster::RasterData& host_mat_c , hs::raster::RasterData& host_mat_d) { /* 输入数据分配显存并拷贝数据 */ hscuda::raster::RasterData device_mat_a; hscuda::raster::RasterData device_mat_b; hscuda::raster::host2device(host_mat_a.view(), device_mat_a); hscuda::raster::host2device(host_mat_b.view(), device_mat_b); /* 输出数据分配显存 */ int size_x = device_mat_a.view().sizeX(); int size_y = device_mat_a.view().sizeY(); hscuda::raster::RasterData device_mat_c(size_x, size_y); hscuda::raster::RasterData device_mat_d(size_x, size_y); /* 执行并行 */ int N = device_mat_a.view().numel(); hs::cuda::elemwise::mimo::Op2in2out(N , device_mat_a.view().pData() , device_mat_b.view().pData() , OddEvenFunctor() , device_mat_c.view().pData() , device_mat_d.view().pData()); cudaDeviceSynchronize(); /* 将结果保存到主机端 */ device_mat_c.toHost(host_mat_c); device_mat_d.toHost(host_mat_d); } ``` --- ### 输入端/输出端逐元素 inputwise/outputwise #### 应用场景 输入端逐元素是个什么意思呢？ 举个例子，**图像的直方图统计**，输出直方图bin个数固定为256，而输出图形的大小不固定，即此时输入和输出数据长度不一致。对于输入数组，需要遍历每一个元素，而对于输出数组，一个元素可能被多次访问。从这个角度，可以看成输入数组的元素要逐次的进行某种独立的运算，因此，就称为“输入端逐元素”。 类似，输出端逐元素的应用场景在**卷积或者相关运算**。 卷积运算输出结果的每一个元素，对应输入数据的多个索引位置的元素，即：对于输入数组，一个元素可能参与多次计算，而对于输出数组，每一个元素的计算又是相互独立的。所以，可以看成输出数组的元素要逐次的进行某种独立的运算，因此，就称为“输出端逐元素”。 进一步推广，element wise 可以看成即满足“输入端逐元素”又满足“输出端逐元素“的特殊情况。 #### 特点 1. 数组长度：输入数组与输出数组长度可以不一致，输入数组与输出数组之间的数据类型没有限制 2. IO配对关系：输入数组与输出数组之间的元素不是一一对应的关系，但是至少在某一端（输出端或输入端）会遍历所有元素逐次进行某种独立运算 3. 在某些情况下，还需要数组的索引作为输入参数。 #### 使用实例 输入端逐元素的模板函数的定义在头文件 [cudalib/public/hs_cuda_inputwise.cuh](cudalib/public/hs_cuda_inputwise.cuh)，相关实例的代码在 [test/cudaInputwise](test/cudaInputwise/) 。 输出端逐元素的模板函数的定义在头文件 [cudalib/public/hs_cuda_outputwise.cuh](cudalib/public/hs_cuda_outputwise.cuh)，相关实例的代码在 [test/cudaOutputwise](test/cudaOutputwise/) 。 ##### 实例(1) 直方图统计 1. 在 [test/cudaInputwise/fun.cu](test/cudaInputwise/fun.cu) 中定义功能算子 ```cpp #include "hs_cuda_inputwise.cuh" struct HistStatisFunctor { HistStatisFunctor(int bin, float* p_hist) : m_bin(bin), m_p_hist(p_hist) { } __device__ void operator()(float val) { int rval = (int)rintf(val); // 四舍五入得到索引 if (rval < m_bin) { // 原子加 atomicAdd(&m_p_hist[rval], 1); } } private: int m_bin; float* m_p_hist; // 输出直返图, 长度为m_bin }; ``` > 说明：我们的输出指针现在是功能算子的成员，因此重载`()`函数不再需要返回值；同时我们也需要定义构造函数来初始化功能算子的各个成员。然后和之前一样，重载`()`函数的形参个数与你要处理的并行计算最小单元个数有关。需要强调的是，成员中的指针必须指向GPU端的地址！ 2. 因为这里只有一个输入（看重载`()`函数），所以使用 `hs::cuda::inputwise::Unary` 模板函数，然后再做进一步的封装： ```cpp void histStatis(const hs::raster::RasterData & mat, int bin, hs::raster::RasterData & hist) { /* 将数据拷贝到的GPU */ hscuda::raster::RasterData device_mat; hscuda::raster::host2device(mat.view(), device_mat); /* 输出分配内存 */ hscuda::raster::RasterData device_hist(bin, 1); /* 调用 inputwise */ int N = device_mat.view().numel(); hs::cuda::inputwise::Unary(N , device_mat.view().pData() , HistStatisFunctor(bin, device_hist.view().pData())); cudaDeviceSynchronize(); /* 输出复制到主机端 */ device_hist.toHost(hist); } ``` ##### 实例(2) padding一圈 0 本实例实现的功能是给定一个矩阵数据，在矩阵的四周填一圈 0，类似这样： ``` img: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 dst_img: 0 0 0 0 0 0 0 0 1 2 3 4 5 0 0 6 7 8 9 10 0 0 11 12 13 14 15 0 0 16 17 18 19 20 0 0 0 0 0 0 0 0 ``` 在这个实例中需要输入数组的索引信息作为功能算子的输入。 1. 在 [test/cudaInputwise/fun.cu](test/cudaInputwise/fun.cu) 中定义功能算子 ```cpp struct PaddingZeroFunctor { PaddingZeroFunctor(int dst_size_x, int dst_size_y , float* p_dst_data , int src_size_x, int src_size_y) : m_dst_size_x(dst_size_x), m_dst_size_y(dst_size_y) , m_p_dst_data(p_dst_data) , m_src_size_x(src_size_x), m_src_size_y(src_size_y) { } __device__ void operator()(int id, float val) { if (id < m_src_size_x * m_src_size_y) { // id 是输入图像的线性索引, 转换为平面索引 int yid = id / m_src_size_x; int xid = id - yid * m_src_size_x; // 加上偏移量 yid = yid + 1; xid = xid + 1; // 转换为输出图线性索引 int dst_id = xid + yid * m_dst_size_x; m_p_dst_data[dst_id] = val; } } private: int m_src_size_y; int m_src_size_x; int m_dst_size_x; int m_dst_size_y; float* m_p_dst_data; }; ``` > 说明：因为需要索引信息，重载`()`函数的第一个形参限定为`int`来表示输入元素的索引值，其余和实例(1)一样的。 2. 因为这里只有一个输入（看重载`()`函数，形参`id`不算），所以使用 `hs::cuda::inputwiseIndex::Unary` 模板函数，然后再做进一步的封装： ```cpp void paddingZero(const hs::raster::RasterData & src_mat, hs::raster::RasterData & dst_mat) { /* 将数据拷贝到GPU */ hscuda::raster::RasterData device_src_mat; hscuda::raster::host2device(src_mat.view(), device_src_mat); /* 输出数据分配内存 */ int size_x = src_mat.view().sizeX(); int size_y = src_mat.view().sizeY(); hscuda::raster::RasterData device_dst_mat(size_x + 2, size_y + 2); // 初始化为0 cudaMemset(device_dst_mat.view().pData(), 0, sizeof(float) * device_dst_mat.view().numel()); /* 调用 inputwiseIndex */ int N = device_src_mat.view().numel(); hs::cuda::inputwiseIndex::Unary(N , device_src_mat.view().pData() , PaddingZeroFunctor( device_dst_mat.view().sizeX() , device_dst_mat.view().sizeY() , device_dst_mat.view().pData() , size_x , size_y) ); cudaDeviceSynchronize(); /* 输出复制到主机端 */ device_dst_mat.toHost(dst_mat); } ``` ##### 实例(3) 栅格数据赋值为固定值 本实例实现功能比较简单，就是将栅格数据的所有元素赋值为一个指定的值。 1. 在 [test/cudaOutputwise/fun.cu](test/cudaOutputwise/fun.cu) 中定义功能算子 ```cpp struct FillValueFunctor { FillValueFunctor(float fill_value) : m_fill_value(fill_value) { } __device__ void operator()(float & out_val) { out_val = m_fill_value; } private: float m_fill_value; }; ``` > 说明：和输入端逐元素不同，现在重载`()`函数的形参得是**引用**，因为形参是输出了，我们需要修改它的值。类似的，输入数据指针现在成为了功能算子的成员，我们需要定义构造函数来初始化各个成员，同时注意数据指针必须指向GPU端的地址。 2. 因为这里只有一个输出（看重载`()`函数），所以使用 `hs::cuda::outputIndex::Unary` 模板函数，然后再做进一步的封装： ```cpp void rasterFillValue(hs::raster::RasterData & mat, float fill_value) { /* 将输入数据拷贝到GPU */ hscuda::raster::RasterData device_mat; hscuda::raster::host2device(mat.view(), device_mat); /* 调用 outputwise 模板函数 */ int N = mat.view().numel(); hs::cuda::outputwise::Unary(N , FillValueFunctor(fill_value) , device_mat.view().pData()); cudaDeviceSynchronize(); /* 将数据拷贝回主机设备 */ device_mat.toHost(mat); } ``` ##### 实例(4) 获取子栅格数据 本实例实现功能是提取栅格数据的某个指定子块，例如： ``` img: 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.00 11.00 12.00 13.00 14.00 15.00 16.00 17.00 18.00 19.00 20.00 21.00 22.00 23.00 24.00 25.00 26.00 27.00 28.00 29.00 30.00 31.00 32.00 33.00 34.00 35.00 36.00 37.00 38.00 39.00 40.00 41.00 42.00 sub img: 9.00 10.00 11.00 16.00 17.00 18.00 23.00 24.00 25.00 30.00 31.00 32.00 ``` 这里的功能算子需要输出数组的索引作为参数 1. 在 [test/cudaOutputwise/fun.cu](test/cudaOutputwise/fun.cu) 中定义功能算子 ```cpp struct SubRasterFunctor { SubRasterFunctor(int offset_x, int offset_y, int sub_size_x, int sub_size_y, int src_size_x, int src_size_y, float* p_src_datas) : m_offset_x(offset_x), m_offset_y(offset_y), m_sub_size_x(sub_size_x), m_sub_size_y(sub_size_y), m_src_size_x(src_size_x), m_src_size_y(src_size_y), m_p_src_datas(p_src_datas) { } __device__ void operator()(int id, float & out_val) { if (id < m_sub_size_x * m_sub_size_y) { // 将输出数据的线性索引转换为平面索引 int sub_yid = id / m_sub_size_x; int sub_xid = id - sub_yid * m_sub_size_x; sub_yid += m_offset_y; sub_xid += m_offset_x; // 将平面索引转换为输入原图像的线性索引 int src_ind = sub_xid + sub_yid * m_src_size_x; if (src_ind < m_src_size_x * m_src_size_y) { out_val = m_p_src_datas[src_ind]; } else { out_val = 0; } } } private: int m_offset_x; int m_offset_y; int m_sub_size_x; int m_sub_size_y; int m_src_size_x; int m_src_size_y; float* m_p_src_datas; }; ``` > 说明：因为需要索引信息，重载`()`函数的第一个形参限定为`int`来表示输入元素的索引值，其余和实例(3)一样的。 2. 因为这里只有一个输出（看重载`()`函数，形参`id`不算），所以使用 `hs::cuda::OutputwiseIndex::Unary` 模板函数，然后再做进一步的封装： ```cpp void getSubRaster(const hs::raster::RasterData & src_mat, int offset_x, int offset_y, int sub_size_x, int sub_size_y, hs::raster::RasterData & dst_mat) { // 将数据拷贝到GPU端 hscuda::raster::RasterData device_src_mat; hscuda::raster::host2device(src_mat.view(), device_src_mat); // 为输出分配显存 hscuda::raster::RasterData device_dst_mat(sub_size_x, sub_size_y); // 调用 outputwiseIndex int src_size_x = src_mat.view().sizeX(); int src_size_y = src_mat.view().sizeY(); int N = device_dst_mat.view().numel(); hs::cuda::outputwiseIndex::Unary(N , SubRasterFunctor(offset_x, offset_y , sub_size_x, sub_size_y , src_size_x, src_size_y , device_src_mat.view().pData()) , device_dst_mat.view().pData() ); cudaDeviceSynchronize(); // 将结果拷贝到主机端 device_dst_mat.toHost(dst_mat); } ``` --- ### 规约操作 Reduction #### 应用场景 规约操作，简单的理解为：给定 N 个输入数据，指定**符合结合律的二元操作符**，最终得到一个规约结果，这个二元操作符便包括求和、求最大值、求最小值、点积等。 #### 特点 1. 数组长度：输入数组的个数限定为 1 个，输出返回一个规约结果值 2. 规约算子的要求：必须为符合结合律的二元操作符，同时需要提供运算初值 3. 有时需要返回对应规约值的索引（如返回最大值所在的索引） #### 使用实例 规约模板函数的定义在头文件 [cudalib/public/hs_cuda_reduce.cuh](cudalib/public/hs_cuda_reduce.cuh)，本实例的代码在 [test/cudaReduce](test/cudaReduce/)。 ##### 实例(1) 规约求和 1. 在 [test/cudaReduce/fun.cu](test/cudaReduce/fun.cu) 中 `#include "hs_cuda_reduce.cuh"` 2. 定义规约二元运算符 ```cpp #include "hs_cuda_reduce.cuh" template struct ReduceSumFunctor { __device__ __forceinline__ T operator()(T a, T b) const { return a + b; } }; ``` 规约算子的要求: 1. 算子的运算必须是符合结合律的二元操作，即重载`()`函数必须为有两个相同类型的形参，并返回当前运算的规约结果 2. 记得使用 const 修饰 3. 调用 `hs_cuda_reduction.cuh` 中的规约模板函数 `hs::cuda::reduce::apply` ```cpp float reduceSum(const hs::raster::RasterData& src_datas) { hscuda::raster::RasterData device_datas; hscuda::raster::host2device(src_datas.view(), device_datas); int N = device_datas.view().numel(); float res_sum = 0.0f; hs::cuda::reduce::apply(N // 输入数据长度 , device_datas.view().pData() // 数据数据指针(设备端) , ReduceSumFunctor() // 二元规约算子 , res_sum // 初值 , &res_sum); // 规约操作结果(返回数据在主机端) cudaDeviceSynchronize(); return res_sum; } ``` ##### 实例(2) 求数组的最大值所在的索引 注意，在本实例中，除了在每次规约运算中保留输出元素，还需要保留对应的索引信息。 1. 在 [test/cudaReduce/fun.cu](test/cudaReduce/fun.cu) 中定义包含**索引信息**的二元规约运算符 ```cpp template struct ReduceMaxIndexFunctor { __device__ __forceinline__ void operator()( T a, int aid // 输入元素a以及对应的索引aid , T b, int bid // 输入元素b以及对应的索引bid , T* output // 当前规约结果 , int* output_index // 当前规约索引 ) const { if (a > b) { (*output) = a; (*output_index) = aid; } else { (*output) = b; (*output_index) = bid; } } }; ``` 2. 调用 `hs_cuda_reduction.cuh` 中带索引的规约模板函数 `hs::cuda::reduceIndex::apply`，并做进一步的封装： ```cpp void reduceMaxIndex(const hs::raster::RasterData& src_datas, float * p_max_val, int * p_max_id) { /* 将输入数据拷贝到的GPU */ hscuda::raster::RasterData device_datas; hscuda::raster::host2device(src_datas.view(), device_datas); int N = device_datas.view().numel(); // 设置规约初值 float init_val = std::numeric_limits::min(); hs::cuda::reduceIndex::apply(N // 输入数据长度 , device_datas.view().pData() // GPU端上的数据指针 , ReduceMaxIndexFunctor() // 带索引的二元规约算子 , init_val // 规约初值 , p_max_val // 输出, 最终的最大值 , p_max_id // 输出, 最终的最大值对应的索引 ); cudaDeviceSynchronize(); } ``` --- ### 扫描操作 Scan #### 应用场景 关于扫描操作，可以参考这篇博客：https://blog.csdn.net/qq_44161734/article/details/134620943 。这里的“扫描”指的就是“前缀和”。 扫描可以分为闭扫描和开扫描两种方式，例如给定一个数组 $a[0], \cdots , a[N-1]$， 闭扫描的计算公式为: $$ c[i] = \sum_{k=0}^i a[k] = a[0] + a[1] + \cdots + a[i]\ ,\ i \in [0, N-1] $$ 开扫描的计算公式为: $$ \begin{aligned} o[0] &= 0\\ o[i] &= \sum_{k=0}^{i-1} a[k] = a[0] + a[1] + \cdots + a[i-1]\ ,\ i \in [1, N-1] \end{aligned} $$ 显然，闭扫描 $c[i]$ 和开扫描 $o[i]$ 之间的关系为: $$ c[i] = o[i] + a[i] $$ > 并行扫描算法在很多算法中都有起到非常重要的基础模块作用，而求前缀和有着非常特殊的作用，详见后面的实例。 #### 特点 在 [cudalib/public/hs_cuda_scan.cuh](cudalib/public/hs_cuda_scan.cuh) 实现的是 [Hillis Steele 算法](https://blog.csdn.net/qq_44161734/article/details/134620943)，最终提供两个模版函数分别实现： ```cpp // 开扫描前缀和 hs::cuda::scan::PrefixOpenSum(int N, const T* p_in, T* p_out); // 闭扫描前缀和 hs::cuda::scan::PrefixSum(int N, const T* p_in, T* p_out); ``` 不同之前的应用场景，scan 不需要自定义功能算子哈，直接用就完了。 #### 使用实例 扫描模板函数的定义在头文件 [cudalib/public/hs_cuda_scan.cuh](cudalib/public/hs_cuda_scan.cuh)，本实例的代码在 [test/cudaScan](test/cudaScan/)。 ##### 实例(1) 求输入数组的开扫描前缀和 比较粗暴，直接上代码 [test/cudaScan/fun.cu](test/cudaScan/fun.cu)，核心就是调用 `hs::cuda::scan::PrefixOpenSum` ```cpp #include "hs_cuda_scan.cuh" void getPrefixOpenSum(const hs::raster::RasterData & raster, hs::raster::RasterData & prefix_sum) { /* 将数据拷贝到GPU端 */ hscuda::raster::RasterData device_raster; hscuda::raster::host2device(raster.view(), device_raster); int N = raster.view().numel(); int size_x = raster.view().sizeX(); int size_y = raster.view().sizeY(); /* 为输出数据分配显存 */ hscuda::raster::RasterData device_prefix_sum(size_x, size_y); /* 直接开扫描调用模板函数即可 */ hs::cuda::scan::PrefixOpenSum(N , device_raster.view().pData() , device_prefix_sum.view().pData()); /* 将结果拷贝到主机端 */ device_prefix_sum.toHost(prefix_sum); } ``` ##### 实例(2) 提取大于5的元素 本实例实现：给定一个矩阵，返回一个其中大于5的所有元素。 这个实例比较综合，需要使用到之前各个实例哈，具体代码见 [test/cudaScan/fun.cu](test/cudaScan/fun.cu)，具体的思路如下： 1. 拷贝输入数据到主机 ```cpp void extractElements( const hs::raster::RasterData & raster // 输入矩阵 , float thre // 阈值, 本实例中等于 5 , hs::raster::RasterData & elems // 输出大于thre的元素 ){ /* 将主机数据复制到GPU */ hscuda::raster::RasterData device_raster; hscuda::raster::host2device(raster.view(), device_raster); // .... ``` 2. 应用 elemwise, 提取大于thre的数据掩膜 ```cpp /* 应用 elemwise, 提取大于thre的数据掩膜 */ hscuda::raster::RasterData mask( device_raster.view().sizeX() , device_raster.view().sizeY()); int N = device_raster.view().numel(); hs::cuda::elemwise::Unary(N , device_raster.view().pData() , GenMaskFunctor(thre) , mask.view().pData() ); cudaDeviceSynchronize(); ``` 功能算子 `GenMaskFunctor` 的实现如下: ```cpp struct GenMaskFunctor { GenMaskFunctor(float thre) : m_thre(thre) { } __device__ float operator()(float val) const { if (val > m_thre) { return 1.0f; } else { return 0.0f; } } private: float m_thre; }; ``` 3. 应用 reduction, 计算掩膜的和, 此值为输出结果的长度 ```cpp float res_sum = 0.0f; hs::cuda::reduce::apply(N , mask.view().pData() , MaskSumFunctor() , res_sum , &res_sum); cudaDeviceSynchronize(); ``` 功能算子 `MaskSumFunctor` 的定义如下： ```cpp struct MaskSumFunctor { __device__ __forceinline__ float operator()(float x, float y) const { return x + y; } }; ``` 4. 应用 scan, 求解掩膜的开扫描前缀和, 开扫描前缀和即是输出元素的索引 ```cpp hscuda::raster::RasterData prefix_sum( mask.view().sizeX() , mask.view().sizeY()); hs::cuda::scan::PrefixOpenSum(N , mask.view().pData() , prefix_sum.view().pData()); cudaDeviceSynchronize(); ``` 5. 应用 inputwise, 遍历输入掩膜和前缀和, 获取掩膜元素为1以及在同样索引位置的前缀和(该前缀和恰好为该元素在输出数组的索引) ```cpp // 输出数据分配显存 hscuda::raster::RasterData device_elems((int)res_sum, 1); hs::cuda::inputwise::Ternary(N , mask.view().pData() , prefix_sum.view().pData() , device_raster.view().pData() , MaskIndexFunctor(device_elems.view().numel() , device_elems.view().pData()) ); cudaDeviceSynchronize(); /* 将结果拷贝到CPU */ device_elems.toHost(elems); ``` 其中，功能算子 `MaskIndexFunctor` 的实现如下： ```cpp struct MaskIndexFunctor { MaskIndexFunctor(int out_size, float* p_out_datas) : m_out_size(out_size), m_p_out_datas(p_out_datas) { } __device__ void operator()(float mask, float fid, float value) { int id = (int)floor(fid); if (id < m_out_size && id >= 0 && mask != 0) { m_p_out_datas[id] = value; } } private: int m_out_size; float* m_p_out_datas; }; ```