【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;
	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

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

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

AI-算力集群通往AGI

背景&#xff1a; 自GPT-4发布以来&#xff0c;全球AI能力的发展势头有放缓的迹象。 但这并不意味着Scaling Law失效&#xff0c;也不是因为训练数据不够&#xff0c;而是结结实实的遇到了算力瓶颈。 具体来说&#xff0c;GPT-4的训练算力约2e25 FLOP&#xff0c;近期发布的几个…

双曲方程初值问题的差分逼近(迎风格式)

稳定性: 数值例子 例一 例二 代码 % function chap4_hyperbolic_1st0rder_1D % test the upwind scheme for 1D hyperbolic equation % u_t + a*u_x = 0,0<x<L,O<t<T, % u(x,0) = |x-1|,0<X<L, % u(0,t) = 1% foundate = 2015-4-22’; % chgedate = 202…

刷代码随想录有感(124):动态规划——最长公共子序列

题干&#xff1a; 代码&#xff1a; class Solution { public:int findLength(vector<int>& nums1, vector<int>& nums2) {vector<vector<int>>dp(nums1.size() 1, vector<int>(nums2.size() 1, 0));int res 0;for(int i 1; i <…

买华为智驾,晚了肯定要后悔

文 | AUTO芯球 作者 | 雷慢 晚了就来不及了&#xff01; 你买华为系的车&#xff0c;薅羊毛真的要趁早。 华为ADS2.0高阶智驾正在慢慢恢复原价&#xff0c; 你看啊&#xff0c;就在昨天&#xff0c;华为宣布ADS智驾优惠后价格调到3万元&#xff0c; 只有6000元的优惠了。…

153. 寻找旋转排序数组中的最小值(中等)

153. 寻找旋转排序数组中的最小值 1. 题目描述2.详细题解3.代码实现3.1 Python3.2 Java 1. 题目描述 题目中转&#xff1a;153. 寻找旋转排序数组中的最小值 2.详细题解 如果不考虑 O ( l o g n ) O(log n) O(logn)的时间复杂度&#xff0c;直接 O ( n ) O(n) O(n)时间复杂…

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

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

vue3实现echarts——小demo

版本&#xff1a; 效果&#xff1a; 代码&#xff1a; <template><div class"middle-box"><div class"box-title">检验排名TOP10</div><div class"box-echart" id"chart1" :loading"loading1"&…

【Portswigger 学院】路径遍历

路径遍历&#xff08;Path traversal&#xff09;又称目录遍历&#xff08;Directory traversal&#xff09;&#xff0c;允许攻击者通过应用程序读取或写入服务器上的任意文件&#xff0c;例如读取应用程序源代码和数据、凭证和操作系统文件&#xff0c;或写入应用程序所访问或…

API Object设计模式

API测试面临的问题 API测试由于编写简单&#xff0c;以及较高的稳定性&#xff0c;许多公司都以不同工具和框架维护API自动化测试。我们基于seldom框架也积累了几千条自动化用例。 •简单的用例 import seldomclass TestRequest(seldom.TestCase):def test_post_method(self…

GDB 远程调试简介

文章目录 1. 前言2. GDB 远程调试2.1 准备工作2.1.1 准备 客户端 gdb 程序2.1.2 准备 服务端 gdbserver2.1.3 准备 被调试程序 2.2 调试2.2.1 通过网络远程调试2.2.1.1 通过 gdbserver 直接启动程序调试2.2.1.2 通过 gdbserver 挂接到已运行程序调试 2.2.2 通过串口远程调试2.2…

紫鸟浏览器搭配IPXProxy代理IP的高效使用指南

​紫鸟指纹浏览器一款专门为跨境电商而生的防关联浏览器&#xff0c;能够帮助跨境电商卖家解决多店铺管理问题。紫鸟指纹浏览器为跨境电商卖家提供稳定的登录环境&#xff0c;并且搭配IP代理&#xff0c;能够解决浏览器指纹记录问题&#xff0c;提高操作的安全性。那如何利用紫…

广州AI绘图模型训练外包定制公司

&#x1f680;设计公司如何借助AI人工智能降本增效&#xff0c;广州这家AI公司值得借鉴— 触站AI&#xff0c;智能图像的创新引擎 &#x1f31f; &#x1f3a8; 触站AI&#xff0c;绘制设计界的未来蓝图 &#x1f3a8;在AI技术的浪潮中&#xff0c;触站AI以其前沿的AI图像技术…

RK3568驱动指南|第十六篇 SPI-第188章 mcp2515驱动编写:复位函数

瑞芯微RK3568芯片是一款定位中高端的通用型SOC&#xff0c;采用22nm制程工艺&#xff0c;搭载一颗四核Cortex-A55处理器和Mali G52 2EE 图形处理器。RK3568 支持4K 解码和 1080P 编码&#xff0c;支持SATA/PCIE/USB3.0 外围接口。RK3568内置独立NPU&#xff0c;可用于轻量级人工…

Redux 使用及基本原理

什么是Redux Redux 是用于js应用的状态管理库&#xff0c;通常和React一起用。帮助开发者管理应用中各个组件之间的状态&#xff0c;使得状态的变化变得更加可预测和易于调试。 Redu也可以不和React组合使用。&#xff08;通常一起使用&#xff09; Redux 三大原则 单一数据源…

在uni-app使用vue3使用vuex

在uni-app使用vue3使用vuex 1.在项目目录中新建一个store目录&#xff0c;并且新建一个index.js文件 import { createStore } from vuex;export default createStore({//数据&#xff0c;相当于datastate: {count:1,list: [{name: 测试1, value: test1},{name: 测试2, value: …

从hugging face 下模型

支持国内下载hugging face 的东西 下模型权重 model_id 是红色圈复制的 代码 记得设置下载的存储位置 import os from pathlib import Path from huggingface_hub import hf_hub_download from huggingface_hub import snapshot_downloadmodel_id"llava-hf/llava-v1…

Swift 中强大的 Key Paths(键路径)机制趣谈(下)

概览 在上一篇博文 Swift 中强大的 Key Paths(键路径)机制趣谈(上)中,我们介绍了 Swift 语言中键路径机制的基础知识,并举了若干例子讨论了它的一些用武之地。 而在本文中我们将再接再厉,继续有趣的键路径大冒险,为 KeyPaths 画上一个圆满的句号。 在本篇博文中,您将…

C++:二维数组的遍历

方式一&#xff1a; #include <vector> #include <iostream> int main() { // 初始化一个2x3的二维向量&#xff08;矩阵&#xff09; std::vector<std::vector<float>> matrix { {1.0, 2.0, 3.0}, // 第一行 {4.0, 5.0, 6.0} // 第二行 };…

企业备份NAS存储一体机

企业文件服务器上的数据、员工电脑里的数据以及NAS存储内数据&#xff0c;需要及时备份&#xff0c;Inforternd存储设备内置了强大的备份服务器功能&#xff0c;无需额外费用&#xff0c;就能轻松将重要数据备份至安全可靠的存储空间中。 无论是GS或GSe 统一存储产品&#xff0…

开放式耳机怎么选?五大2024年口碑销量爆棚机型力荐!

在选购开放式耳机的时候&#xff0c;我们总会因为有太多的选择而陷入两难&#xff0c;又想要一个颜值比较高的&#xff0c;又想要同时兼顾性能还不错的&#xff0c;所以作为测评博主&#xff0c;今天我们就给大家带来自己的一些选购技巧和自己觉得还不错开放式耳机&#xff0c;…