【CUDA】 归约 Reduction

news2025/1/12 3:57:02

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;
	else
		std::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;
	else
		std::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;
	else
		std::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;
	else
		std::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;
	else
		std::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 memory
    partial_sum[tx] = i < n ? X[i] : 0;
    __syncthreads();

    // Accumulate the elements using reduction
    for(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 memory
    partial_sum[tx] = i < n ? X[i] : 0;
    __syncthreads();

    // Accumulate the elements using reduction
    for(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 memory
    partial_sum[tx] = i < n ? X[i] : 0;
    __syncthreads();

    // Accumulate the elements using reduction
    for(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 memory
    partial_sum[tx] = i < n ? X[i] : 0;
    partial_sum[tx] += i + blockDim.x < n ? X[i + blockDim.x] : 0;
    __syncthreads();

    // Accumulate the elements using reduction
    for(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 memory
    partial_sum[tx] = i < n ? X[i] : 0;
    partial_sum[tx] += i + blockDim.x < n ? X[i + blockDim.x] : 0;
    __syncthreads();    

    // Accumulate the elements using reduction
    for(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.coloradmin.cn/o/1888973.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

三菱A系列网络连接

寄存器名 读写 寄存器类型 变量类型 寄存器范围 说明 X##1 R/W BIT I/O离散 0&#xff0d;7FF Input Y##1 R/W BIT I/O离散 0&#xff0d;7FF Output M##1 R/W BIT I/O离散 0&#xff0d;9255 Internal relay B##1 R/W BIT I/O离散 0&#xff0d;3FF Link relay F##1 R/W BIT I…

EPS绘制甘家寨地形图

1、数据准备 &#xff08;1&#xff09;外业采集的数据点&#xff1b; &#xff08;2&#xff09;地形草图 2、软件准备 这里准备的是EPS2021版本的绘图软件&#xff0c;如下&#xff1a; 3、开始绘图 &#xff08;1&#xff09;打开软件&#xff0c;如上图&#xff0c;选择【…

不同行业如何选择适合自己行业的项目管理工具?

在当今的信息化时代&#xff0c;项目管理软件已成为各行各业不可或缺的工具。然而&#xff0c;由于各行业具有不同的特点和需求&#xff0c;因此选择合适的项目管理软件成为了一个重要问题。本文将探讨不同行业在选择项目管理软件时需要考虑的因素&#xff0c;希望能帮助大家更…

8.12 矢量图层面要素单一符号使用十四(标记符号渲染边界)

前言 本章介绍矢量图层线要素单一符号中标记符号渲染边界&#xff08;Outline: Marker line&#xff09;的使用说明&#xff1a;文章中的示例代码均来自开源项目qgis_cpp_api_apps 标记符号渲染边界&#xff08;Outline: Marker line&#xff09; Outline系列只画边界&#…

HarmonyOS ArkUi 官网踩坑:单独隐藏导航条无效

环境&#xff1a; 手机&#xff1a;Mate 60 Next版本&#xff1a; NEXT.0.0.26 导航条介绍 导航条官网设计指南 setSpecificSystemBarEnabled 设置实际效果&#xff1a; navigationIndicator&#xff1a;隐藏导航条无效status&#xff1a;会把导航条和状态栏都隐藏 官方…

深入理解策略梯度算法

策略梯度&#xff08;Policy Gradient&#xff09;算法是强化学习中的一种重要方法&#xff0c;通过优化策略以获得最大回报。本文将详细介绍策略梯度算法的基本原理&#xff0c;推导其数学公式&#xff0c;并提供具体的例子来指导其实现。 策略梯度算法的基本概念 在强化学习…

Ajax异步请求 axios

Ajax异步请求 axios 1 axios介绍 原生ajax请求的代码编写太过繁琐,我们可以使用axios这个库来简化操作&#xff01; 在后续学习的Vue(前端框架)中发送异步请求,使用的就是axios. 需要注意的是axios不是vue的插件,它可以独立使用. axios说明网站&#xff1a;(https://www.kancl…

猫头虎分享[可灵AI」官方推荐的驯服指南-V1.0

猫头虎分享[可灵AI」官方推荐的驯服指南-V1.0 猫头虎是谁&#xff1f; 大家好&#xff0c;我是 猫头虎&#xff0c;别名猫头虎博主&#xff0c;擅长的技术领域包括云原生、前端、后端、运维和AI。我的博客主要分享技术教程、bug解决思路、开发工具教程、前沿科技资讯、产品评…

云服务器安装部署LAMP网站Web环境教程

搭建网站如何安装LAMP环境&#xff0c;以腾讯云轻量应用服务器为例&#xff0c;应用模板直接选择“LAMP”镜像即可&#xff0c;打开腾讯云轻量应用服务器页面&#xff0c;在应用模板中选择LAMP即可&#xff0c;如下图&#xff1a; 轻量服务器“LAMP”镜像 腾讯云的LAMP应用镜像…

XLSX + LuckySheet + LuckyExcel实现前端的excel预览

文章目录 功能简介简单代码实现效果参考 功能简介 通过LuckyExcel的transformExcelToLucky方法&#xff0c; 我们可以把一个文件直接转成LuckySheet需要的json字符串&#xff0c; 之后我们就可以用LuckySheet预览excelLuckyExcel只能解析xlsx格式的excel文件&#xff0c;因此对…

海南云亿商务咨询有限公司抖音电商新引擎

在当今这个数字化高速发展的时代&#xff0c;电商已经成为推动经济增长的重要引擎。而在众多电商平台中&#xff0c;抖音以其独特的短视频形式和庞大的用户群体&#xff0c;迅速崛起为电商领域的新星。海南云亿商务咨询有限公司&#xff0c;作为一家专注于抖音电商服务的公司&a…

JavaScript实现注册页面的校验

完成注册页面的校验 1.目标 要求&#xff1a; 1、 用户名必须是6-10位的字母或者数字 2、 密码长度必须6位任意字符 3、 两次输入密码要一致 说明&#xff1a;只要有一个输入项不满足要求则阻止表单提交。都满足才可以提交表单。 2.实现 1.知识点 1.1js事件 【1】鼠标离…

网上最靠谱的改名大师颜廷利:世界顶级哲学家,东方伟大的思想家教育家

颜廷利教授&#xff0c;21世纪的东方哲学巨匠、科学探索者&#xff0c;以及中国教育界的领军人物&#xff0c;他在《升命学说》中阐述了一种深刻的生活哲学。这本书包含了四个主要理念&#xff1a;净化论、和合法则、唯悟主义与镜正理念&#xff0c;每一个都为我们如何理解生活…

我在手提电脑上将大模型训练成了语文老师

&#xff08;图片由大模型生成&#xff0c;如有侵权&#xff0c;立删&#xff09; 记得一年多以前&#xff0c;和不少商家交流大模型解决方案时&#xff0c;他们谈到内部有很多的资料&#xff0c;可以对大模型进行训练&#xff0c;让大模型变得更有智慧&#xff0c;从而为客户…

【项目日记(三)】搜索引擎-搜索模块

❣博主主页: 33的博客❣ ▶️文章专栏分类:项目日记◀️ &#x1f69a;我的代码仓库: 33的代码仓库&#x1f69a; &#x1faf5;&#x1faf5;&#x1faf5;关注我带你了解更多项目内容 目录 1.前言2.项目回顾3.搜索流程3.1分词3.2触发3.3去重3.4排序3.5包装 4.总结 1.前言 在前…

全网新鲜出炉的Stable Diffusion 人物发型提示词大全,中英文列表!

前言 简介&#xff1a; 使用发型提示词能更精确描述所需图像的发型特征&#xff0c;如卷发、短发、颜色和风格。结合正负提示词&#xff0c;确保生成图片符合预期。尝试使用工具如PromptChoose来创建个性化图像描述&#xff0c;包含多种发型选项&#xff0c;如刘海、马尾、波浪…

6.5、函数的常见形式

代码 #include <iostream> using namespace std; #include <string>//函数的的常见延时 //1、无参无反 void test01() {cout << "this is test01" << endl; } //2、有参无反 void test02(int a) {cout << "this is test02 a &q…

QT学习积累——方法参数加const和不加const的区别

目录 引出方法参数加const和不加const的区别方法加static和不加static的区别Qt遍历list提高效率显示函数的调用使用&与不使用&除法的一个坑 总结自定义信号和槽1.自定义信号2.自定义槽3.建立连接4.进行触发 自定义信号重载带参数的按钮触发信号触发信号拓展 lambda表达…

【探索Linux】P.36(传输层 —— TCP协议段格式)

阅读导航 引言一、TCP段的基本格式二、控制位详细介绍三、16位接收窗口大小⭕窗口大小的作用⭕窗口大小的限制⭕窗口缩放选项⭕窗口大小的更新⭕窗口大小与拥塞控制 四、紧急指针温馨提示 引言 在上一篇文章中&#xff0c;我们深入探讨了一种无连接的UDP协议&#xff0c;它以其…

【数据结构】04.双向链表

一、双向链表的结构 注意&#xff1a;这里的“带头”跟前面我们说的“头节点”是两个概念&#xff0c;带头链表里的头节点&#xff0c;实际为“哨兵位”&#xff0c;哨兵位节点不存储任何有效元素&#xff0c;只是站在这里“放哨的”。 “哨兵位”存在的意义&#xff1a;遍历循…