在上一篇文章中，我们介绍了如何使用GPU运行的并行算法。这些并行任务是相互之间完全独立的任务，这与我们通常所知道的编程方式有很大的不同。虽然我们可以从并行中获益，但这种奇妙的并行操作方式一定让我们非常感兴趣。复杂的。所以在本文的Numba代码中，我们将介绍一些让线程协同计算的常用技术。首先加载相应的包fromtimeimportperf_counterimportnumpyasnpimportnumbafromnumbaimportcudaprint(np.__version__)print(numba.__version__)cuda.detect()#1.21.6#0.55.2#找到1个CUDA设备#id0b'TeslaT4'[支持]#计算能力：7.5#PCI设备ID：4#PCI总线ID：0#UUID：GPU-bcc6196e-f26e-afdc-1db3-6eba6ff160f0#看门狗：禁用#FP32/FP64PerformanceRatio:32#Summary:#1/1devicesaresupported#True别忘了，我们是来CUDA编程的，所以NV的GPU是必须的，比如colab或者Kaggle可以免费上。线程之间的协调简单的并行归约算法我们将从一个非常简单的问题开始本节：对数组的所有元素求和。该算法非常简单。如果不使用NumPy，我们可以这样实现它：defsum_cpu(array):s=0.0foriinrange(array.size):s+=array[i]returns这看起来不太像Pythonic。但它确实让我们知道它正在跟踪数组中的所有元素。如果s的结果取决于数组的每个元素，我们如何并行化该算法？首先，我们需要重写算法以允许并行化，如果有部分不能并行化，应该允许线程相互通信。到目前为止，我们还没有学会如何让线程之间进行通信……其实我们之前说过，不同块中的线程是不通信的。我们可以考虑只启动一个块，但正如我们上次所说，一个块在大多数GPU上只能有1024个线程！如何克服这个？如果将数组拆分为1024个块（或适当数量的threads_per_block）并分别对每个块求和怎么办？然后在最后，我们可以对每个块求和的结果进行求和。下图显示了一个非常简单的2块拆分示例。上图是对数组元素求和的“分而治之”法。这如何在GPU上完成？首先，您需要将数组拆分成块。每个数组块将恰好对应一个具有固定线程数的CUDA块。在每个块中，每个线程可以对多个数组元素求和。然后将这些per-thread值相加，需要线程进行通信，我们将在下一个例子中讨论。由于我们正在并行化块，因此内核的输出应该设置为一个块。为了减少，我们将它复制到CPU并在那里完成工作。threads_per_block=1024#为什么不呢！blocks_per_grid=32*80#使用32*倍数的流式多处理器#示例2.1：原始缩减@cuda.jitdefreduce_naive(array,partial_reduction):i_start=cuda.grid(1)threads_per_grid=cuda.blockDim.x*cuda.gridDim.xs_thread=0.0fori_arrinrange(i_start,array.size,threads_per_grid):s_thread+=array[i_arr]#我们需要创建一个特殊的*共享*数组，它可以被读取和写入块中的每个线程。每个块都有自己的共享数组。请参阅下面的警告！s_block=cuda.shared.array((threads_per_block,),numba.float32)#我们现在将单个线程的本地临时总和存储到共享数组中。由于共享数组的大小为#threads_per_block==blockDim.x#（本例中为1024），我们应该使用`threadIdx.x`对其进行索引。tid=cuda.threadIdx.xs_block[tid]=s_thread#下一行同步块中的线程。它确保在#该行之后，所有值都已写入`s_block`。cuda.syncthreads()#最后，我们需要对所有线程的值求和，以产生一个#每个块的值。为此，我们只需要一个线程。iftid==0:#我们将共享数组元素的总和存储在它的第一个#坐标中foriinrange(1,threads_per_block):s_block[0]+=s_block[i]#将这个部分总和移动到输出。这里只有一个线程在写。partial_reduction[cuda.blockIdx.x]=s_block[0]这里需要注意的是数组必须是共享的，每个数组块都取“small”这里的“small”：具体大小取决于数组的计算能力GPU，通常在48KB和163KB之间请参阅此表项中的“每个线程块的最大共享内存量”。在编译时有一个已知大小（这就是我们调整共享数组threads_per_block而不是blockDim.x的原因）。我们总是可以为任何大小的共享数组定义一个工厂函数……但要注意这些内核的编译时间。这里的数组需要为Numba类型指定的数据类型，而不是Numpy类型（没有理由！）。N=1_000_000_000a=np.arange(N,dtype=np.float32)a/=a.sum()#a的sum=1（达到float32精度）s_cpu=a.sum()#高度优化的NumPyCPU代码timing_cpu=np.empty(21)foriinrange(timing_cpu.size):tic=perf_counter()a.sum()toc=perf_counter()timing_cpu[i]=toc-tictiming_cpu*=1e3#转换为msprint(f"运行时间CPU:{timing_cpu.mean():.0f}±{timing_cpu.std():.0f}ms")#运行时间CPU:354±24msdev_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()#CPU的最终缩减np.isclose(s,s_cpu)#确保我们有正确的数字#Truetiming_naive=np.empty(21)foriinrange(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()assertnp.isclose(s,s_cpu)timing_naive[i]=toc-tictiming_naive*=1e3#转换为msprint(f"Elapsedtimenaive:{timing_naive.mean():.0f}±{timing_naive.std():.0f}ms")#Elapsedtimenaive:30±12ms在GoogleColab上测试，得到了10倍的加速题外话：上面的方法被称为简单归约算法，因为它是最简单和最容易实现的。我们在大数据中常见的Map-Reduce算法就是这个算法。虽然实现简单，但易于理解，因此很常用。当然，它也是出了名的慢。有兴趣的可以研究一下。更好的并行缩减算法上面的算法是最“幼稚”的，所以有很多技巧可以加速这种代码（参见CUDA演示文稿中的优化并行缩减以获取基准）。在介绍更好的方法之前，让我们回顾一下内核函数的最后一点：iftid==0:#Singlethreadtakingcareofbusinessforiinrange(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()。#示例2.2：更好的缩减@cuda.jitdefreduce_better(array,partial_reduction):i_start=cuda.grid(1)threads_per_grid=cuda.blockDim.x*cuda.gridDim.xs_thread=0.0fori_arrinrange(i_start,array.size,threads_per_grid):s_thread+=array[i_arr]#我们需要创建一个特殊的*共享*数组，块中的每个线程都可以读取和写入该数组。每个块都有自己的共享数组。请参阅下面的警告！s_block=cuda.shared.array((threads_per_block,),numba.float32)#我们现在将线程的本地临时总和存储到共享数组中。#由于共享数组的大小为threads_per_block==blockDim.x，#我们应该使用`threadIdx.x`对其进行索引。tid=cuda.threadIdx.xs_block[tid]=s_thread#下一行同步块中的线程。它确保在#该行之后，所有值都已写入s_block。cuda同步器hreads()i=cuda.blockDim.x//2while(i>0):if(tid
