从头开始进行CUDA编程：线程间协作的常见技术

转载：Deephub Imba

在前一篇文章中，我们介绍了如何使用 GPU 运行的并行算法。这些并行任务是那些完全相互独立的任务，这点与我们一般认识的编程方式有很大的不同，虽然我们可以从并行中受益，但是这种奇葩的并行运行方式对于我们来说肯定感到非常的复杂。所以在本篇文章的Numba代码中，我们将介绍一些允许线程在计算中协作的常见技术。

首先还是载入相应的包

from time import perf_counterimport numpy as npimport numbafrom numba import cudaprint(np.__version__)print(numba.__version__)cuda.detect()# 1.21.6# 0.55.2# Found 1 CUDA devices# id 0             b'Tesla T4'                              [SUPPORTED]#                       Compute Capability: 7.5#                            PCI Device ID: 4#                               PCI Bus ID: 0#                                     UUID: GPU-bcc6196e-f26e-afdc-1db3-6eba6ff160f0#                                 Watchdog: Disabled#              FP32/FP64 Performance Ratio: 32# Summary:# 1/1 devices are supported# True

不要忘记，我们这里是CUDA编程，所以NV的GPU是必须的，比如可以去colab或者Kaggle白嫖。

线程间的协作

简单的并行归约算法

我们将从一个非常简单的问题开始本节:对数组的所有元素求和。这个算法非常简单。如果不使用NumPy，我们可以这样实现它:

def sum_cpu(array):s = 0.0for i in range(array.size):s += array[i]return s

这看起来不是很 Pythonic。但它能够让我们了解它正在跟踪数组中的所有元素。如果 s 的结果依赖于数组的每个元素，我们如何并行化这个算法呢？首先，我们需要重写算法以允许并行化，如果有无法并行化的部分则应该允许线程相互通信。

到目前为止，我们还没有学会如何让线程相互通信……事实上，我们之前说过不同块中的线程不通信。我们可以考虑只启动一个块，但是我们上次也说了，在大多数 GPU 中块只能有 1024 个线程！

如何克服这一点？如果将数组拆分为 1024 个块（或适当数量的threads_per_block）并分别对每个块求和呢？然后最后，我们可以将每个块的总和的结果相加。下图显示了一个非常简单的 2 块拆分示例。

上图就是对数组元素求和的“分而治之”方法。

如何在 GPU 上做到这一点呢？首先需要将数组拆分为块。每个数组块将只对应一个具有固定数量的线程的CUDA块。在每个块中，每个线程可以对多个数组元素求和。然后将这些每个线程的值求和，这里就需要线程进行通信，我们将在下一个示例中讨论如何通信。

由于我们正在对块进行并行化，因此内核的输出应该被设置为一个块。为了完成Reduce，我们将其复制到 CPU 并在那里完成工作。

threads_per_block = 1024  # Why not!blocks_per_grid = 32 * 80  # Use 32 * multiple of streaming multiprocessors# Example 2.1: Naive reduction@cuda.jitdef reduce_naive(array, partial_reduction):i_start = cuda.grid(1)threads_per_grid = cuda.blockDim.x * cuda.gridDim.xs_thread = 0.0for i_arr in range(i_start, array.size, threads_per_grid):s_thread += array[i_arr]# We need to create a special *shared* array which will be able to be read# from and written to by every thread in the block. Each block will have its# own shared array. See the warning below!s_block = cuda.shared.array((threads_per_block,), numba.float32)# We now store the local temporary sum of a single the thread into the# shared array. Since the shared array is sized#     threads_per_block == blockDim.x# (1024 in this example), we should index it with `threadIdx.x`.tid = cuda.threadIdx.xs_block[tid] = s_thread# The next line synchronizes the threads in a block. It ensures that after# that line, all values have been written to `s_block`.cuda.syncthreads()# Finally, we need to sum the values from all threads to yield a single# value per block. We only need one thread for this.if tid == 0:# We store the sum of the elements of the shared array in its first# coordinatefor i in range(1, threads_per_block):s_block[0] += s_block[i]# Move this partial sum to the output. Only one thread is writing here.partial_reduction[cuda.blockIdx.x] = s_block[0]

这里需要注意的是必须共享数组，并且让每个数组块变得“小”

这里的“小”：确切大小取决于 GPU 的计算能力，通常在 48 KB 和 163 KB 之间。请参阅此表中的“每个线程块的最大共享内存量”项。（https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#features-and-technical-specifications__technical-specifications-per-compute-capability）

在编译时有一个已知的大小（这就是我们调整共享数组threads_per_block而不是blockDim.x的原因）。我们总是可以为任何大小的共享数组定义一个工厂函数……但要注意这些内核的编译时间。

这里的数组需要为 Numba 类型指定的 dtype，而不是 Numpy 类型（这个没有为什么！）。

N = 1_000_000_000a = np.arange(N, dtype=np.float32)a /= a.sum() # a will have sum = 1 (to float32 precision)s_cpu = a.sum()# Highly-optimized NumPy CPU codetiming_cpu = np.empty(21)for i in range(timing_cpu.size):tic = perf_counter()a.sum()toc = perf_counter()timing_cpu[i] = toc - tictiming_cpu *= 1e3  # convert to msprint(f"Elapsed time CPU: {timing_cpu.mean():.0f} ± {timing_cpu.std():.0f} ms")# Elapsed time CPU: 354 ± 24 msdev_a = cuda.to_device(a)dev_partial_reduction = cuda.device_array((blocks_per_grid,), dtype=a.dtype)reduce_naive[blocks_per_grid, threads_per_block](dev_a, dev_partial_reduction)s = dev_partial_reduction.copy_to_host().sum()  # Final reduction in CPUnp.isclose(s, s_cpu)  # Ensure we have the right number# Truetiming_naive = np.empty(21)for i in range(timing_naive.size):tic = perf_counter()reduce_naive[blocks_per_grid, threads_per_block](dev_a, dev_partial_reduction)s = dev_partial_reduction.copy_to_host().sum()cuda.synchronize()toc = perf_counter()assert np.isclose(s, s_cpu)    timing_naive[i] = toc - tictiming_naive *= 1e3  # convert to msprint(f"Elapsed time naive: {timing_naive.mean():.0f} ± {timing_naive.std():.0f} ms")# Elapsed time naive: 30 ± 12 ms

在谷歌Colab上测试，得到了10倍的加速。

题外话：上面这个方法之所以说是简单的规约算法，是因为这个算法最简单，也最容易实现。我们在大数据中常见的Map-Reduce算法就是这个算法。虽然实现简单，但是他容易理解，所以十分常见，当然他慢也是出名的，有兴趣的大家可以去研究研究。

一种更好的并行归约算法

上面的算法最 “朴素”的，所以有很多技巧可以加快这种代码的速度（请参阅 CUDA 演示文稿中的 Optimizing Parallel Reduction 以获得基准测试）。

在介绍更好的方法之前，让我们回顾一下内核函数的的最后一点：

if tid == 0:  # Single thread taking care of businessfor i in range(1, threads_per_block):s_block[0] += s_block[i]partial_reduction[cuda.blockIdx.x] = s_block[0]

我们并行化了几乎所有的操作，但是在内核的最后，让一个线程负责对共享数组 s_block 的所有 threads_per_block 元素求和。为什么不能把这个总和也并行化呢？

听起来不错对吧，下图显示了如何在 threads_per_block 大小为 16 的情况下实现这一点。我们从 8 个线程开始工作，第一个将对 s_block[0] 和 s_block[8] 中的值求和。第二个求和s_block[1]和s_block[9]中的值，直到最后一个线程将s_block[7]和s_block[15]的值相加。

在下一步中，只有前 4 个线程需要工作。第一个线程将对 s_block[0] 和 s_block[4] 求和；第二个，s_block[1] 和 s_block[5]；第三个，s_block[2] 和 s_block[6]；第四个也是最后一个，s_block[3] 和 s_block[7]。

第三步，只需要 2 个线程来处理 s_block 的前 4 个元素。

第四步也是最后一步将使用一个线程对 2 个元素求和。

由于工作已在线程之间分配，因此它是并行化的。这不是每个线程的平等划分，但它是一种改进。在计算上，这个算法是 O(log2(threads_per_block))，而上面“朴素”算法是 O(threads_per_block)。比如在我们这个示例中是 1024 次操作，用于了两个算法差距有10倍

最后还有一个细节。在每一步，我们都需要确保所有线程都已写入共享数组。所以我们必须调用 cuda.syncthreads()。

# Example 2.2: Better reduction@cuda.jitdef reduce_better(array, partial_reduction):i_start = cuda.grid(1)threads_per_grid = cuda.blockDim.x * cuda.gridDim.xs_thread = 0.0for i_arr in range(i_start, array.size, threads_per_grid):s_thread += array[i_arr]# We need to create a special *shared* array which will be able to be read# from and written to by every thread in the block. Each block will have its# own shared array. See the warning below!s_block = cuda.shared.array((threads_per_block,), numba.float32)# We now store the local temporary sum of the thread into the shared array.# Since the shared array is sized threads_per_block == blockDim.x,# we should index it with `threadIdx.x`.tid = cuda.threadIdx.xs_block[tid] = s_thread# The next line synchronizes the threads in a block. It ensures that after# that line, all values have been written to `s_block`.cuda.syncthreads()i = cuda.blockDim.x // 2while (i > 0):if (tid < i):s_block[tid] += s_block[tid + i]cuda.syncthreads()i //= 2if tid == 0:partial_reduction[cuda.blockIdx.x] = s_block[0]reduce_better[blocks_per_grid, threads_per_block](dev_a, dev_partial_reduction)s = dev_partial_reduction.copy_to_host().sum()  # Final reduction in CPUnp.isclose(s, s_cpu)# Truetiming_naive = np.empty(21)for i in range(timing_naive.size):tic = perf_counter()reduce_better[blocks_per_grid, threads_per_block](dev_a, dev_partial_reduction)s = dev_partial_reduction.copy_to_host().sum()cuda.synchronize()toc = perf_counter()assert np.isclose(s, s_cpu)    timing_naive[i] = toc - tictiming_naive *= 1e3  # convert to msprint(f"Elapsed time better: {timing_naive.mean():.0f} ± {timing_naive.std():.0f} ms")# Elapsed time better: 22 ± 1 ms

可以看到，这比原始方法快25%。

重要说明：你可能很想将同步线程移动到 if 块内，因为在每一步之后，超过当前线程数一半的内核将不会被使用。但是这样做会使调用同步线程的 CUDA 线程停止并等待所有其他线程，而所有其他线程将继续运行。因此停止的线程将永远等待永远不会停止同步的线程。如果您同步线程，请确保在所有线程中调用 cuda.syncthreads()。

i = cuda.blockDim.x // 2while (i > 0):if (tid < i):s_block[tid] += s_block[tid + i]#cuda.syncthreads()  这个是错的cuda.syncthreads()  # 这个是对的i //= 2

Numba 自动归约

其实归约算法并不简单，所以Numba 提供了一个方便的 cuda.reduce 装饰器，可以将二进制函数转换为归约。所以上面冗长而复杂的算法可以替换为：

@cuda.reducedef reduce_numba(a, b):return a + b# Compile and checks = reduce_numba(dev_a)np.isclose(s, s_cpu)# True# Timetiming_numba = np.empty(21)for i in range(timing_numba.size):tic = perf_counter()s = reduce_numba(dev_a)toc = perf_counter()assert np.isclose(s, s_cpu)    timing_numba[i] = toc - tictiming_numba *= 1e3  # convert to msprint(f"Elapsed time better: {timing_numba.mean():.0f} ± {timing_numba.std():.0f} ms")# Elapsed time better: 45 ± 0 ms

上面的运行结果我们可以看到手写代码通常要快得多（至少 2 倍），但 Numba 给我们提供的方法却非常容易使用。这对我们来说是格好事，因为终于有我们自己实现的Python方法比官方的要快了

从头开始进行CUDA编程：线程间协作的常见技术相关推荐

Java并发编程—线程间协作方式wait()、notify()、notifyAll()和Condition
原文作者:Matrix海子原文地址:Java并发编程:线程间协作的两种方式:wait.notify.notifyAll和Condition 目录一.wait().notify()和notifyA ...
LinuxC高级编程——线程间同步
LinuxC高级编程--线程间同步宗旨:技术的学习是有限的,分享的精神是无限的. 1. 互斥锁mutex 多个线程同时访问共享数据时可能会冲突.对于多线程的程序,访问冲突的问题是很普遍的,解决的办法 ...
19、Java并发编程：线程间协作的两种方式：wait、notify、notifyAll和Condition
Java并发编程:线程间协作的两种方式:wait.notify.notifyAll和Condition 在前面我们将了很多关于同步的问题,然而在现实中,需要线程之间的协作.比如说最经典的生产者-消费者 ...
线程间协作的两种方式：wait、notify、notifyAll和Condition
转载自线程间协作的两种方式:wait.notify.notifyAll和Condition 在前面我们将了很多关于同步的问题,然而在现实中,需要线程之间的协作.比如说最经典的生产者-消费者模型:当 ...
JAVA线程间协作：wait.notify.notifyAll
欢迎支持笔者新作:<深入理解Kafka:核心设计与实践原理>和<RabbitMQ实战指南>,同时欢迎关注笔者的微信公众号:朱小厮的博客. 欢迎跳转到本文的原文链接:https: ...
Linux多线程编程---线程间同步（互斥锁、条件变量、信号量和读写锁）
本篇博文转自http://zhangxiaoya.github.io/2015/05/15/multi-thread-of-c-program-language-on-linux/ Linux下提供了 ...
Java并发——线程间通信与同步技术
传统的线程间通信与同步技术为Object上的wait().notify().notifyAll()等方法,Java在显示锁上增加了Condition对象,该对象也可以实现线程间通信与同步.本文会介绍有 ...
C++并发编程线程间共享数据std::future和sd::promise
线程间共享数据使用互斥锁实现线程间共享数据为了避免死锁可以考虑std::lock()或者boost::shared_mutex 要尽量保护更少的数据同步并发操作 C++标准库提供了一些工具可以 ...
《Java-SE-第二十四章》之线程间协作
前言在你立足处深挖下去,就会有泉水涌出!别管蒙昧者们叫嚷:"下边永远是地狱!" 博客主页:KC老衲爱尼姑的博客主页博主的github,平常所写代码皆在于此共勉:talk is ...

从头开始进行CUDA编程：线程间协作的常见技术

线程间的协作

从头开始进行CUDA编程：线程间协作的常见技术相关推荐

最新文章

热门文章