【CUDA】 归约 Reduction

Reduction

Reduction算法从一组数值中产生单个数值。这个单个数值可以是所有元素中的总和、最大值、最小值等。

图1展示了一个求和Reduction的例子。

图1

线程层次结构

在Reduction算法中，线程的常见组织方式是为每个元素使用一个线程。下面将展示利用许多不同方式来组织线程。

Code

Host代码用随机值初始化输入向量，并调用kernel执行Reduction。

#include <iostream>
#include <cstdio>
#include <ctime>
#include <cmath>
#include <cuda_runtime.h>
#include <thrust/host_vector.h>
#include <thrust/device_vector.h>//float array[6] = { 3, 1, 2, 3, 5, 4 };
//float* dev_array = 0;
//cudaMalloc(&dev_array, 4 * 6);
//cudaMemcpy(dev_array, array, 4 * 6, cudaMemcpyHostToDevice);
//thrust::device_ptr<float> dev_ptr(dev_array);
//thrust::reduce(dev_ptr, dev_ptr + 6);//由于dev_array指向device端，不能直接作为参数，需要对其封装
//
//thrust::host_vector<type> hvec;
//thrust::device_vector<type> dvec;
//dvec = hvec;
//thrust::reduce(dvec.begin(), dvec.end());//此时的参数是迭代器，不用也不能用device_ptr对其封装
//
上述的两种函数的调用方法也存在host端的版本，传入的指针或者迭代器都是host端数据
//thrust::reduce(array, array + 6);
//thrust::reduce(hvec.begin(), hvec.end());
//
从device_ptr中提取“原始”指针需要使用raw_pointer_cast函数
//float dev_array = thrust::raw_pointer_cast(dev_ptr);#include "helper_cuda.h"
#include "error.cuh"
#include "reductionKernels.cuh"const double EPSILON = 1.0;
const int TIMENUMS = 50;int main(void)
{int n, t_x;n = 4096;t_x = 1024;thrust::host_vector<double> h_A(n);for (int i = 0; i < n; ++i) {h_A[i] = rand() / (double)RAND_MAX;}double h_res = thrust::reduce(h_A.begin(), h_A.end(), 0.0f, thrust::plus<double>());int temp = n;thrust::device_vector<double> d_A;//device vector 和 host vector 可以直接用等号进行传递，对应于cudaMemcpy的功能thrust::device_vector<double> d_B = h_A;dim3 threads(t_x);dim3 gridDim;cudaEvent_t start, stop;float elapsed_time;float sumTime = 0;for (int i = 0; i < TIMENUMS; i++) {temp = n;d_B = h_A;do {gridDim = dim3((temp + threads.x - 1) / threads.x);d_A = d_B;d_B.resize(gridDim.x);checkCudaErrors(cudaEventCreate(&start));checkCudaErrors(cudaEventCreate(&stop));checkCudaErrors(cudaEventRecord(start));reduction1<double> << <gridDim, threads, threads.x * sizeof(double) >> > (thrust::raw_pointer_cast(d_A.data()),temp,thrust::raw_pointer_cast(d_B.data()));checkCudaErrors(cudaEventRecord(stop));checkCudaErrors(cudaEventSynchronize(stop));checkCudaErrors(cudaEventElapsedTime(&elapsed_time, start, stop));sumTime += elapsed_time;temp = gridDim.x;//std::cout << "\t" << h_res << " != " << d_B[0] << "\t" << d_B[1]  << std::endl;} while (temp > 1);}printf("Time = %g ms.\n", sumTime / TIMENUMS);if (abs(h_res - d_B[0]) > 0.01)std::cout << "Error (Reduction 1):" << h_res << " != " << d_B[0] << std::endl;elsestd::cout << "Reduction 1: Success!" << std::endl;sumTime = 0;for (int i = 0; i < TIMENUMS; i++) {temp = n;d_B = h_A;do {gridDim = dim3((temp + threads.x - 1) / threads.x);d_A = d_B;d_B.resize(gridDim.x);checkCudaErrors(cudaEventRecord(start));reduction2<double> << <gridDim, threads, threads.x * sizeof(double) >> > (thrust::raw_pointer_cast(d_A.data()),temp,thrust::raw_pointer_cast(d_B.data()));checkCudaErrors(cudaEventRecord(stop));checkCudaErrors(cudaEventSynchronize(stop));checkCudaErrors(cudaEventElapsedTime(&elapsed_time, start, stop));sumTime += elapsed_time;temp = gridDim.x;} while (temp > 1);}printf("Time = %g ms.\n", sumTime / TIMENUMS);if (abs(h_res - d_B[0] > 0.01))std::cout << "Error (Reduction 2): " << h_res << " != " << d_B[0] << std::endl;elsestd::cout << "Reduction 2: Success!" << std::endl;sumTime = 0;for (int i = 0; i < TIMENUMS; i++) {temp = n;d_B = h_A;do {gridDim = dim3((temp + threads.x - 1) / threads.x);d_A = d_B;d_B.resize(gridDim.x);checkCudaErrors(cudaEventRecord(start));reduction3<double> << <gridDim, threads, threads.x * sizeof(double) >> > (thrust::raw_pointer_cast(d_A.data()),temp,thrust::raw_pointer_cast(d_B.data()));checkCudaErrors(cudaEventRecord(stop));checkCudaErrors(cudaEventSynchronize(stop));checkCudaErrors(cudaEventElapsedTime(&elapsed_time, start, stop));sumTime += elapsed_time;temp = gridDim.x;} while (temp > 1);}printf("Time = %g ms.\n", sumTime / TIMENUMS);if (abs(h_res - d_B[0] > 0.01))std::cout << "Error (Reduction 3): " << h_res << " != " << d_B[0] << std::endl;elsestd::cout << "Reduction 3: Success!" << std::endl;sumTime = 0;for (int i = 0; i < TIMENUMS; i++) {temp = n;d_B = h_A;//checkCudaErrors(cudaEventRecord(start));do {gridDim = dim3((temp + threads.x - 1) / threads.x);d_A = d_B;d_B.resize(gridDim.x);checkCudaErrors(cudaEventRecord(start));reduction4<double> << <gridDim, threads, threads.x * sizeof(double) >> > (thrust::raw_pointer_cast(d_A.data()),temp,thrust::raw_pointer_cast(d_B.data()));checkCudaErrors(cudaEventRecord(stop));checkCudaErrors(cudaEventSynchronize(stop));checkCudaErrors(cudaEventElapsedTime(&elapsed_time, start, stop));sumTime += elapsed_time;temp = gridDim.x;} while (temp > 1);}printf("Time = %g ms.\n", sumTime / TIMENUMS);if (abs(h_res - d_B[0] > 0.01))std::cout << "Error (Reduction 4): " << h_res << " != " << d_B[0] << std::endl;elsestd::cout << "Reduction 4: Success!" << std::endl;sumTime = 0;for (int i = 0; i < TIMENUMS; i++) {temp = n;d_B = h_A;do {gridDim = dim3((temp + threads.x - 1) / threads.x);d_A = d_B;d_B.resize(gridDim.x);checkCudaErrors(cudaEventRecord(start));reduction5<double> << <gridDim, threads, threads.x * sizeof(double) >> > (thrust::raw_pointer_cast(d_A.data()),temp,thrust::raw_pointer_cast(d_B.data()));checkCudaErrors(cudaEventRecord(stop));checkCudaErrors(cudaEventSynchronize(stop));checkCudaErrors(cudaEventElapsedTime(&elapsed_time, start, stop));sumTime += elapsed_time;temp = gridDim.x;} while (temp > 1);}printf("Time = %g ms.\n", sumTime / TIMENUMS);if (abs(h_res - d_B[0] > 0.01))std::cout << "Error (Reduction 5): " << h_res << " != " << d_B[0] << std::endl;elsestd::cout << "Reduction 5: Success!" << std::endl;return 0;
}

Note：

helper_cuda.h 与error.cuh头文件为错误查验工具。后续会发布到Github。

去除checkCudaErrors等错误查验函数不影响程序运行。

Reduction 1: Interleaved Addressing - Diverge Warps

template<typename T> __global__
void reduction1(T* X, uint32_t n, T* Y){extern __shared__ uint8_t shared_mem[];T* partial_sum = reinterpret_cast<T*>(shared_mem);uint32_t tx = threadIdx.x;uint32_t i = blockIdx.x * blockDim.x + threadIdx.x;// Load the elements of the array into the shared memorypartial_sum[tx] = i < n ? X[i] : 0;__syncthreads();// Accumulate the elements using reductionfor(uint32_t stride = 1; stride < blockDim.x; stride <<= 1){if(tx % (2 * stride) == 0)partial_sum[tx] += tx + stride < n ? partial_sum[tx + stride] : 0;__syncthreads();}if(tx == 0) Y[blockIdx.x] = partial_sum[0];
}

首先，kernel将数组的元素加载到共享内存中。然后，在每次迭代中，它将上次迭代的相邻元素相加。第一次迭代会添加相邻的元素。第二次迭代会添加相隔2个元素的元素，依此类推。当步幅大于块中线程数时，循环结束。

最后，kernel将结果写入输出数组。如图2所示：

图2

使用这个kernel，添加元素的线程不是连续的，而且在每次迭代中，线程之间的差距会更大。这将导致warp需要重新运行大量的代码。

Reduction 2: Interleaved Addressing - Bank Conflicts

template<typename T> __global__
void reduction2(T* X, uint32_t n, T* Y){extern __shared__ uint8_t shared_mem[];T* partial_sum = reinterpret_cast<T*>(shared_mem);uint32_t tx = threadIdx.x;uint32_t i = blockIdx.x * blockDim.x + threadIdx.x;// Load the elements of the array into the shared memorypartial_sum[tx] = i < n ? X[i] : 0;__syncthreads();// Accumulate the elements using reductionfor(uint32_t stride = 1; stride < blockDim.x; stride <<= 1){int index = 2 * stride * tx;if (index < blockDim.x)partial_sum[index] += index + stride < n ? partial_sum[index + stride] : 0;__syncthreads();}if(tx == 0) Y[blockIdx.x] = partial_sum[0];
}

这个kernel和以前那个很相似。唯一的区别是添加元素的线程是连续的。这样可以让线程更有效地运行代码。如图3

图3

这个kernel的问题在于它引发了Bank Conflicts（板块冲突）。当两个线程访问共享内存的同一个存储区时就会发生Bank Conflicts。计算能力在2.0及以上的设备有32个32位宽的存储区。如果两个线程访问相同的存储区，第二次访问将延迟直到第一次访问完成。

Reduction 3: Sequential Addressing

template<typename T> __global__
void reduction3(T* X, uint32_t n, T* Y){extern __shared__ uint8_t shared_mem[];T* partial_sum = reinterpret_cast<T*>(shared_mem);uint32_t tx = threadIdx.x;uint32_t i = blockIdx.x * blockDim.x + threadIdx.x;// Load the elements of the array into the shared memorypartial_sum[tx] = i < n ? X[i] : 0;__syncthreads();// Accumulate the elements using reductionfor(uint32_t stride = blockDim.x / 2; stride > 0; stride >>= 1){if(tx < stride)partial_sum[tx] += partial_sum[tx + stride];__syncthreads();}if(tx == 0) Y[blockIdx.x] = partial_sum[0];
}

为了避免Bank Conflicts，这个kernel使用顺序寻址。在每次迭代中，线程和元素都是连续的。如图4所示：

这个kernel的问题是，一半的线程在第一次循环迭代中是空闲的。

图4

Reduction 4: First Add During Load

template<typename T> __global__
void reduction4(T* X, uint32_t n, T* Y){extern __shared__ uint8_t shared_mem[];T* partial_sum = reinterpret_cast<T*>(shared_mem);uint32_t tx = threadIdx.x;uint32_t i = blockIdx.x * (blockDim.x * 2) + threadIdx.x;// Load the elements of the array into the shared memorypartial_sum[tx] = i < n ? X[i] : 0;partial_sum[tx] += i + blockDim.x < n ? X[i + blockDim.x] : 0;__syncthreads();// Accumulate the elements using reductionfor(uint32_t stride = blockDim.x / 2; stride > 0; stride >>= 1){if(tx < stride)partial_sum[tx] += partial_sum[tx + stride];__syncthreads();}if(tx == 0) Y[blockIdx.x] = partial_sum[0];
}

这个kernel和Reduction 3的kernel类似。唯一的区别是在迭代开始之前，每个线程加载并添加数组的两个元素。

Reduction 5: Unroll Last Warp

template<typename T> __device__ 
void warpReduce(volatile T* partial_sum, uint32_t tx){partial_sum[tx] += partial_sum[tx + 32];partial_sum[tx] += partial_sum[tx + 16];partial_sum[tx] += partial_sum[tx +  8];partial_sum[tx] += partial_sum[tx +  4];partial_sum[tx] += partial_sum[tx +  2];partial_sum[tx] += partial_sum[tx +  1];
}template<typename T, int block_size> __global__
void reduction5(T* X, uint32_t n, T* Y){extern __shared__ uint8_t shared_mem[];T* partial_sum = reinterpret_cast<T*>(shared_mem);uint32_t tx = threadIdx.x;uint32_t i = blockIdx.x * (blockDim.x * 2) + threadIdx.x;// Load the elements of the array into the shared memorypartial_sum[tx] = i < n ? X[i] : 0;partial_sum[tx] += i + blockDim.x < n ? X[i + blockDim.x] : 0;__syncthreads();    // Accumulate the elements using reductionfor(uint32_t stride = blockDim.x / 2; stride > 32; stride >>= 1){if(tx < stride)partial_sum[tx] += partial_sum[tx + stride];__syncthreads();}if(tx < 32) warpReduce(partial_sum, tx);if(tx == 0) Y[blockIdx.x] = partial_sum[0];
}

kernel再次与Reduction 4相似。唯一的区别是最后一个warp是展开的。这将使warp能够在无控制分歧的情况下运行代码。

性能分析

运行时间：

向量维度：4096

线程块维度：1024

由运行时间分析，前三种算法耗时逐次递减。Reduction 4、5，较Reduction 3，虽然再加载共享内存时多计算了一次，但是引入了更多的全局内存访问，所以Reduction 4表现与Reduction 3相比略差，但是同时Reduction 4、5会减少设备与host的传输次数，当传输数据量大时Reduction 4、5整体上可能有更好的表现。Reduction 5相比Reduction 4，减少了线程束的重复执行，所以耗时略有减少。

Note：单次运行可能因为设备启动原因，各种算法运行时间差异较大，可采用循环20次以上取平均值。

笔者采用设备：RTX3060 6GB

PMPP项目提供的分析

kernel的性能是使用NvBench项目在多个gpu中测量的。研究的性能测量方法有:

内存带宽:每秒传输的数据量。

内存带宽利用率:占用内存带宽的百分比。

Reduction 1: Interleaved Addressing - Diverge Warps

Reduction 2: Interleaved Addressing - Bank Conflicts

Reduction 3: Sequential Addressing

Reduction 4: First Add During Load

Reduction 5: Unroll Last Warp

参考文献：

1、大规模并行处理器编程实战（第2版）

2、​​​PPMP

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.mzph.cn/news/865410.shtml

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈email:809451989@qq.com，一经查实，立即删除！

相关文章

反射快速入门

反射快速入门

反射就是通过字节码文件获取类的成员变量、构造方法和成员方法的所有信息。 利用反射&#xff0c;我们可以获取成员变量的修饰符、名字、类型、取值。我们可以获取构造方法的名字、形参&#xff0c;并利用通过反射获取的构造方法创建对象。我们可以获取成员方法的修饰符、名字、…
阅读更多...
主食冻干复查｜希喂、喜崽、生生不息可以盲选吗？测评结果来揭秘

主食冻干复查｜希喂、喜崽、生生不息可以盲选吗？测评结果来揭秘

在挑选主食冻干时&#xff0c;许多宠物主人都会感到头疼。尽管主食冻干相较于普通猫粮具有诸多优势&#xff0c;但其价格也相对高昂。这导致许多宠物主人担心高价购买的主食冻干可能营养价值并不理想。然而&#xff0c;在选择时&#xff0c;我们还需要考虑其他重要因素&#xf…
阅读更多...
Zabbix 配置PING监控

Zabbix 配置PING监控

Zabbix PING监控介绍 如果需要判断机房的网络或者主机是否正常&#xff0c;这就需要使用zabbix ping&#xff0c;Zabbix使用外部命令fping处理ICMP ping的请求&#xff0c;在基于ubuntu APT方式安装zabbix后默认已存在fping程序。另外zabinx_server配置文件参数FpingLocation默…
阅读更多...
【C++】多态(详解)

【C++】多态(详解)

前言&#xff1a;今天学习的内容可能是近段时间最难的一个部分的内容了&#xff0c;C的多态&#xff0c;这部分内容博主认为难度比较大&#xff0c;各位一起慢慢啃下来。 &#x1f496; 博主CSDN主页:卫卫卫的个人主页 &#x1f49e; &#x1f449; 专栏分类:高质量&#xff23…
阅读更多...
操作系统真象还原：编写硬盘驱动程序

操作系统真象还原：编写硬盘驱动程序

第13章-编写硬盘驱动程序 这是一个网站有所有小节的代码实现&#xff0c;同时也包含了Bochs等文件 13.1 硬盘及分区表 13.1.1 创建从盘及获取安装的磁盘数 要实现文件系统&#xff0c;必须先有个磁盘介质&#xff0c;虽然咱们己经有个虚拟磁盘 hd60M.img&#xff0c;但它只…
阅读更多...
力扣67 二进制求和

力扣67 二进制求和

文章目录 1. 题目链接2. 题目代码3.感受 1. 题目链接 二进制求和 2. 题目代码 class Solution { public:string addBinary(string a, string b) {vector<int> stringA;vector<int> stringB;int lengthOfA a.length();int lengthOfB b.length();for(int subscrip…
阅读更多...
OceanBase Meetup北京站｜跨行业应用场景中的一体化分布式数据库：AI赋能下的探索与实践

OceanBase Meetup北京站｜跨行业应用场景中的一体化分布式数据库：AI赋能下的探索与实践

随着业务规模的不断扩张和数据处理需求的日益复杂化&#xff0c;传统数据库架构逐渐暴露出业务稳定性波动、扩展性受限、处理效率降低以及运营成本高等一系列问题。众多行业及其业务场景纷纷踏上了数据库现代化升级之路。 为应对这些挑战&#xff0c;7月6日&#xff0c;OceanB…
阅读更多...
专注于文件夹加密和保护的免费软件

专注于文件夹加密和保护的免费软件

一、简介 1、这是一款专注于文件夹加密和保护的免费软件。允许用户为重要的文件或文件夹设置密码&#xff0c;从而防止未经授权的访问。软件提供了隐藏、锁定、只读等多种保护模式&#xff0c;用户可以根据需要选择适合的模式来保护文件。除了基本的加密功能外&#xff0c;它还…
阅读更多...
【java计算机毕设】陪诊师管理系统java MySQL springboot vue3 Maven源码 代码+文档PPT

【java计算机毕设】陪诊师管理系统java MySQL springboot vue3 Maven源码 代码+文档PPT

目录 1项目功能 2项目介绍 3项目地址 1项目功能 【java计算机毕设】陪诊师管理系统java MySQL springboot vue3 Maven项目设计源码代码万字文档ppt 2项目介绍 系统功能&#xff1a; vue3陪诊师管理系统。 该平台采用了前后端分离技术&#xff0c;SpringBoot和VUE3框架&…
阅读更多...
告别熬夜改稿：AI降重工具让论文降重变得轻松又有趣

告别熬夜改稿：AI降重工具让论文降重变得轻松又有趣

已经天临五年了&#xff0c;大学生们还在为论文降重烦恼……手动降重确实是个难题&#xff0c;必须要先付点小经费去靠谱的网站查重&#xff0c;再对着红字标注去改&#xff0c;后面每一次的论文呢查重结果都像赌//博&#xff0c;谁也不知道明明是同一篇文章&#xff0c;第二次…
阅读更多...
Halcon 曲线追踪

Halcon 曲线追踪

Halcon 曲线追踪&#xff08;边缘检测、xld分割、xld筛选、线段合并&#xff09; 图片数据与程序 链接&#xff1a;https://pan.baidu.com/s/1feGOa0A7dvCeBjQivr6TvA 提取码&#xff1a;f2ws 原图 起点终点方向 * 1.加载图片 ********************************************…
阅读更多...
Python处理异常用操作介绍

Python处理异常用操作介绍

Python中的异常处理主要用于捕获和处理程序运行过程中出现的错误。 在编写Python程序时&#xff0c;我们经常会遇到各种错误&#xff0c;如语法错误、运行时错误等。为了确保程序的稳定性和健壮性&#xff0c;我们需要对可能出现的错误进行捕获和处理。本文将介绍Python中常用的…
阅读更多...
[笔记] 卷积 - 02 滤波器在时域的等效形式

[笔记] 卷积 - 02 滤波器在时域的等效形式

1.讨论 这里主要对时域和频域的卷积运算的特征做了讨论&#xff0c;特别是狄拉克函数的物理意义。 关于狄拉克函数&#xff0c;参考这个帖子&#xff1a;https://zhuanlan.zhihu.com/p/345809392 1.狄拉克函数提到的好函数的基本特征是能够快速衰减&#xff0c;对吧&#xf…
阅读更多...
软件功能测试基础知识大揭秘，功能测试报告就找专业软件测评机构

软件功能测试基础知识大揭秘，功能测试报告就找专业软件测评机构

软件功能测试是以软件产品的需求规格为基础&#xff0c;通过对软件功能的逐个测试&#xff0c;验证软件是否符合需求规格&#xff0c;是否能够正常执行各项功能操作。对于软件产品而言&#xff0c;功能测试是一项至关重要的工作&#xff0c;它能够发现软件中存在的功能缺陷、错…
阅读更多...
多微信运营管理方案

多微信运营管理方案

微信作为一款社交通讯软件&#xff0c;已经成为人们日常生活中不可缺少的工具。不仅个人&#xff0c;很多企业都用微信来联系客户、维护客户和营销&#xff0c;这自然而然就会有很多微信账号、手机也多&#xff0c;那管理起来就会带来很多的不便&#xff0c;而多微信私域管理系…
阅读更多...
softmax从零开始实现

softmax从零开始实现

softmax从零开始实现 代码结果 代码 import numpy as np import torch import torchvision import torchvision.transforms as transforms from torch.utils import data# H,W,C -> C,H,W mnist_train torchvision.datasets.FashionMNIST(root"./data", trainTr…
阅读更多...
java静态代理-被代理对象，代理对象的概念（图+代码解释）

java静态代理-被代理对象，代理对象的概念（图+代码解释）

案例是老师类&#xff0c;这个老师生病请假了&#xff0c;需要请另外一个老师临时帮忙&#xff0c;这个过来帮忙的老师就是代理对象&#xff0c;生病的老师就是被代理对象&#xff0c;其中我们需要代理对象和被代理对象都implement这个ITeacherDao接口&#xff0c;实现里面的te…
阅读更多...
8款你不一定知道的良心软件！

8款你不一定知道的良心软件！

AI视频生成&#xff1a;小说文案智能分镜智能识别角色和场景批量Ai绘图自动配音添加音乐一键合成视频https://aitools.jurilu.com/我们使用一些流行的软件的时候&#xff0c;往往会忽略一些功能非常强大的软件&#xff0c;因为这些软件的众 多&#xff0c;都因为看不见而丢失&a…
阅读更多...
udp发送数据如果超过1个mtu时，抓包所遇到的问题记录说明

udp发送数据如果超过1个mtu时，抓包所遇到的问题记录说明

最近在测试Syslog udp发送相关功能&#xff0c;测试环境是centos udp头部的数据长度是2个字节&#xff0c;最大传输长度理论上是65535&#xff0c;除去头部这些字节&#xff0c;可以大概的说是64k。 写了一个超过64k的数据(随便用了一个7w字节的buffer)发送demo&#xff0c;打…
阅读更多...
从百数教学看产品设计：掌握显隐规则，打造极致用户体验

从百数教学看产品设计：掌握显隐规则，打造极致用户体验

字段显隐规则允许通过一个控件&#xff08;如复选框、单选按钮或下拉菜单&#xff09;来控制其他控件&#xff08;如文本框、日期选择器等&#xff09;和标签页&#xff08;如表单的不同部分&#xff09;的显示或隐藏。 这种规则通常基于用户的选择或满足特定条件来触发&#…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023