使用 Cuda 并行实现计算大型数组中的大型连续子序列之和

Posted 2023-04-15

【中文标题】使用 Cuda 并行实现计算大型数组中的大型连续子序列之和

【英文标题】：Parallel implementation of the computation of the sum of large contiguous subsequences in a large array using Cuda

【发布时间】：2017-07-02 11:11:19

【问题描述】：

我有一个包含 1000 个元素的大型数组，我想使用 CUDA 计算这个大型数组中大小为 100 的大型连续子序列的总和。 这是一个小尺寸的说明性示例。数组大小为20，序列大小为5。

tab = [80,12,14,5,70,9,26,30,8,12,16,15,60,12,38,32,17,67,19,11]

顺序如下：

S1= 80+12+14+5+70
S2= 12+14+5+70+9
S3= 14+5+70+9+26
....

您是否有一个有效的想法来使用具有 1000 个元素的数组和 100 个序列的 CUDA 并行化此任务？

【问题讨论】：

前几天你不是问了几乎完全相同的问题吗？这有什么不同？

事实上，前一个是针对小数组和小序列的，建议的解决方案不适用于具有大序列的大数组。我在上一个问题中提到，在我的情况下，我使用的是一个大数组，但它已经足够清晰了。@talonemies。

1000 个元素不是一个大数组。但是 - 你有多少序列？一百？一千？一百万？如果是后一种情况，那么我们正在谈论。

【参考方案1】：

一些序言：

正如 cmets 中已经提到的，1000 的总数组大小对于在 GPU 上运行是一个非常小问题。如果这是您在 GPU 上做的唯一事情，那么您可能不会从 GPU 代码中获得有趣的加速。 previous answer 建议使用一维模板。即使模板宽度为 100（或 200），仍然应该可以使用 1D 模板方法。这里唯一真正限制模板大小的因素是共享内存的大小，共享内存可以轻松保存等于每个块的线程数 + 模板（序列）“光环”的 200 的数据大小。但是该主题已经在此处介绍，因此我将介绍另一种方法，它应该能够处理 任意 数组长度和 任意 序列大小。

在这种情况下，看看我们是否可以利用prefix sum 来帮助我们可能是合适的。前缀和存在快速并行方法，因此这很有吸引力。前缀和只是一个数字的输出序列，表示输入序列中所有先前数字的总和。所以如果我有一个这样的数组：

1, 1, 3, 2, 0, 1

独占前缀总和将是：

0, 1, 2, 5, 7, 7, 8

排除这里的意思是当前求和位置包含所有之前的值，但排除当前位置的数据项。对于排他前缀和，请注意，我们可以（如果我们愿意）为比输入序列长度长 1 的序列生成“有用”数据。

前缀和可以很容易地适应您所问的问题。假设对于我上面的序列，我们想要计算 3 的子序列。我们可以取前缀和，并从中减去相同的前缀和序列，将 右 移动 3（序列长度），例如所以：

  0, 1, 2, 5, 7, 7, 8
-          0, 1, 2, 5, 7, 7, 8
=          5, 6, 5, 3

在这种情况下，这个序列 (5, 6, 5, 3) 将是所需的答案。

下面是一个完整的例子，只是在推力上。我之所以这样做，是因为在 CUDA 中编写快速并行前缀和并不是一件容易的事，因此我更愿意使用（并建议其他人使用）library 实现。如果您想探索如何编写自己的并行前缀和，可以探索开源的推力，或阅读this 的介绍。

这是一个完整的例子：

$ cat t1279.cu
#include <thrust/device_vector.h>
#include <thrust/scan.h>
#include <thrust/transform.h>
#include <thrust/functional.h>
#include <thrust/copy.h>
#include <iostream>

const int arr_size =  1000;
const int seq_len  =   100;

typedef int mytype;

int main()

  // test case
  mytype t_arr[] = 80,12,14,5,70,9,26,30,8,12,16,15,60,12,38,32,17,67,19,11,0;
  int t_arr_size = sizeof(t_arr)/sizeof(mytype);
  int t_seq_len = 5;
  thrust::device_vector<mytype> d_arr1(t_arr, t_arr+t_arr_size);
  thrust::device_vector<mytype> d_res1(t_arr_size);
  thrust::device_vector<mytype> d_out1(t_arr_size-t_seq_len);

  thrust::exclusive_scan(d_arr1.begin(), d_arr1.end(), d_res1.begin());
  thrust::transform(d_res1.begin()+t_seq_len, d_res1.end(), d_res1.begin(), d_out1.begin(), thrust::minus<mytype>());
  thrust::copy_n(d_out1.begin(), t_arr_size-t_seq_len, std::ostream_iterator<mytype>(std::cout, ","));
  std::cout << std::endl;

  // case with larger array length and larger sequence length
  thrust::device_vector<mytype> d_arr(arr_size+1, 1);
  thrust::device_vector<mytype> d_res(arr_size+1);
  thrust::device_vector<mytype> d_out(arr_size+1-seq_len);

  thrust::inclusive_scan(d_arr.begin(), d_arr.end(), d_res.begin());
  thrust::transform(d_res.begin()+seq_len, d_res.end(), d_res.begin(), d_out.begin(), thrust::minus<mytype>());
  // validate
  for (int i = 0; i < arr_size+1-seq_len; i++) 
    mytype t = d_out[i];
    if (t != seq_len) std::cout << "mismatch at: " << i << "was: " << t << "should be: " << seq_len << std::endl; return 1;
    
  return 0;


$ nvcc -arch=sm_35 -o t1279 t1279.cu
$ ./t1279
181,110,124,140,143,85,92,81,111,115,141,157,159,166,173,146,
$

【讨论】：

所以这具有 O(n) 复杂度，而不是一维模板的 O(n*m)，其中 m 为 100 或模板大小？

非常感谢您的回答，亲爱的罗伯特先生。事实上，我的实习导师让我从这个开始，然后我将转向更大的数据。我仍然想知道是否有可能与前一个循环平行，或者只是我忘记了它？我是 CUDA 和并行性的新手，如果我问了很多问题，我很抱歉。@Robert Crovella

我不确定您的意思是“与前一个循环平行”。我收集到您正试图通过在线程之间重用尽可能多的中间和来最小化计算量。我想说，实际上，一个写得很好的并行前缀和就是这样做的。如果您阅读the material I indicated，您可能会得出相同的结论。

【参考方案2】：

非常相似的问题：Parallel implementation of the computation of the sum of contiguous subsequences in an array using Cuda

罗伯特·克罗维拉回答。

这在序列 size=4 之前得到了回答，现在这只是 size=100 的不同。要么让内核适应这个类似的问题，要么在 cpu 上解决这个问题，因为 cpu 应该能够在 1000-10000 个周期内解决这个（N=1000，L=100）问题（复杂度为 O(N) ）核心，而 gpu 的许多周期将未被使用，只有高级 gpu 可以进行负载平衡以使用这些空闲资源。 CPU 甚至可以在将数组发送到 gpu 之前解决这个问题。但如果这必须在 GPU 上发生（在 2 个内核之间），那么我只能建议对 Robert Crovella 的答案进行与硬件相关的优化，例如：

从纹理缓存和/或常量内存中获取数组的一半，以进一步增加内存带宽（在共享/全局内存之上）提供计算资源。（还必须能够以这种方式增加每个计算单元/smx 的可用大小) 使用多个临时累加器而不是一个，以增加编译器执行某些指令优化的机会（指令级并行性？） 如果 FP64 资源的数量与 FP32 资源的数量相当，则使用双精度（或单精度）来增加更多性能，但仅在 N 类似于 2 * total_FP32_units / total_FP64_units 的每 N 个线程组上，并且仅当有多个每个 smx 的 warps 以确保 FP32 和 FP64 资源都被使用。也许其中一个累加器可以是 FP64，而不是在整个组上使用 FP64。（但这会隐藏更多的转换延迟） 不仅计算第 i 个索引，还计算第 i+1st 和 i+2nd 个索引，每个线程只需要每个线程额外添加 4 次，而不是在计算第 i 次后增加 200 次（在另外两个线程上）一。例如，添加下一个元素并减去最旧的元素。为什么每个线程有 3 个作业？保持所有内存条忙碌，而不是进行内存条冲突（如果总共有 2、4、8、10 个内存条）。如果有 3、6、9 个内存库，那么它可能是 2 个或 4 个作业线程。但这会增加（双倍、三倍）临时内存需求，因此必须有一个性能平衡点，可以对某些架构进行基准测试以了解它。

我不知道cuda，但它应该能够使用多个内存限定符来使用更多的gpu芯片区域来更快地提供数据，就像opencl一样。

如果共享内存和全局内存有不同/单独的管道到达流处理器，那么您甚至可以将全局内存用于某些工作（我的意思是，不将所有内容加载到共享，而是将一些从全局加载，并且让它使用所有 shared_lines + global_lines)

32M 元素数组可能无法从第一次优化中获得帮助，但您仍然可以使用 FP64 累加器来拥有多个临时累加器（例如，每个线程使用 8 个 FP32 累加器和 1 个 FP64 累加器（使用 FMA 来减少其延迟可能是隐藏在这 8 个累加器后面），所以整个数组元素被平等对待）但是如果常量内存或纹理缓存有 50000 个元素的空间，它可以服务 250 个组（每个组有 200 个来自这种类型的内存的元素）或 1250 个组（每个有50 个来自此内存，150 个来自其他内存（共享，..））并且可能仍会稍微提高性能。（例如，100k 元素数组需要 1k 个组，%25 的线程使用常量内存）

【讨论】：

每次迭代都会重用之前的总和？@huseyin tugrul buyukisik

我完全改变了答案，如果是你，你应该重新考虑你的支持。

谢谢@huseyin tugrul buyukisik 的回答，我对这些东西真的很陌生，我会仔细阅读你的回答并努力做到。

@Robert Crovella 罗伯特先生有必要的设备来尝试这个 FP64 东西，也许他会尝试并与我们分享结果。如果他能说出真相，我将不胜感激。

我希望@Robert Crovella