我使用 C++ 和 CUDA/C 并想为特定问题编写代码，但遇到了一个非常棘手的缩减问题。
我在并行编程方面的经验并非微不足道，但非常有限，我不能完全预见这个问题的特殊性。我怀疑是否有一种方便甚至“简单”的方法来处理我面临的问题，但也许我错了。如果有任何资源（即文章、书籍、网络链接……）或涵盖此问题或类似问题的关键字，请告诉我。
我试图尽可能地概括整个案例并保持抽象而不是发布太多代码。
布局 ...
我有一个由 N 个初始元素和 N 个结果元素组成的系统。（例如，我将使用 N=8，但 N 可以是任何大于三的整数值。）
static size_t const N = 8;
double init_values[N], result[N];
我需要计算几乎所有（恐怕不是全部）初始值的唯一排列，而不会产生自干扰。
这意味着计算f(init_values[0],init_values[1]), f(init_values[0],init_values[2]), ..., f(init_values[0],init_values[N-1]), f(init_values[1],init_values[2]), ..., f(init_values[1],init_values[N-1]), ... 等等。
这实际上是一个虚拟三角矩阵，其形状如下图所示。
P 0 1 2 3 4 5 6 7
|---------------------------------------
0| x
|
1| 0 x
|
2| 1 2 x
|
3| 3 4 5 x
|
4| 6 7 8 9 x
|
5| 10 11 12 13 14 x
|
6| 15 16 17 18 19 20 x
|
7| 21 22 23 24 25 26 27 x
每个元素都是 中相应列和行元素的函数init_values。
P[i] (= P[row(i)][col(i]) = f(init_values[col(i)], init_values[row(i)])
IE
P[11] (= P[5][1]) = f(init_values[1], init_values[5])
使用示例有(N*N-N)/2 = 28可能的独特组合（注意：P[1][5]==P[5][1]，所以我们只有一个下（或上）三角矩阵）N = 8。
基本问题
结果数组从 P 计算为行元素的总和减去相应列元素的总和。例如，位置 3 的结果将计算为第 3 行的总和减去第 3 列的总和。
result[3] = (P[3]+P[4]+P[5]) - (P[9]+P[13]+P[18]+P[24])
result[3] = sum_elements_row(3) - sum_elements_column(3)
我试图用 N = 4 的图片来说明它。
因此，以下情况属实：
N-1操作（潜在的并发写入）将在每个result[i]result[i]将从减法和加法中N-(i+1)写入i- 每次输出
P[i][j]都会有一个减法r[j]和一个加法r[i]
这是主要问题出现的地方：
- 使用一个线程计算每个 P 并直接更新结果将导致多个内核尝试写入相同的结果位置（每个 N-1 个线程）。
- 另一方面，存储整个矩阵 P 用于后续的归约步骤在内存消耗方面非常昂贵，因此对于非常大的系统来说是不可能的。
为每个线程块拥有一个独特的、共享的结果向量的想法也是不可能的。（50k 的 N 产生 25 亿个 P 元素，因此 [假设每个块最多有 1024 个线程] 如果每个块都有自己的结果数组和 50k 双元素，则最少 240 万个块消耗超过 900GiB 的内存。）
我认为我可以处理更静态行为的减少，但就潜在的并发内存写访问而言，这个问题是相当动态的。（或者是否可以通过一些“基本”类型的减少来处理它？）
增加一些并发症...
不幸的是，根据与初始值无关的（任意用户）输入，需要跳过 P 的某些元素。假设我们需要跳过排列 P[6]、P[14] 和 P[18]。因此，我们还剩下 24 个组合，需要计算。
如何告诉内核需要跳过哪些值？我提出了三种方法，如果 N 非常大（如数万个元素），每种方法都有明显的缺点。
1.存储所有组合...
...具有各自的行和列索引struct combo { size_t row,col; };，需要在 a 中计算vector<combo>并对该向量进行操作。（由当前实现使用）
std::vector<combo> elements;
// somehow fill
size_t const M = elements.size();
for (size_t i=0; i<M; ++i)
{
// do the necessary computations using elements[i].row and elements[i].col
}
该解决方案消耗大量内存，因为只有“几个”（甚至可能是一万个元素，但与总共数十亿个元素相比并不多），但它避免了
- 索引计算
- 找到移除的元素
对于 P 的每个元素，这是第二种方法的缺点。
2.对P的所有元素进行操作，找到移除的元素
如果我想对 P 的每个元素进行操作并避免嵌套循环（我无法在 cuda 中很好地重现），我需要执行以下操作：
size_t M = (N*N-N)/2;
for (size_t i=0; i<M; ++i)
{
// calculate row indices from `i`
double tmp = sqrt(8.0*double(i+1))/2.0 + 0.5;
double row_d = floor(tmp);
size_t current_row = size_t(row_d);
size_t current_col = size_t(floor(row_d*(ict-row_d)-0.5));
// check whether the current combo of row and col is not to be removed
if (!removes[current_row].exists(current_col))
{
// do the necessary computations using current_row and current_col
}
}
与第一个示例中的向量相比，该向量removes非常小，elements但是获得 和 if 分支的额外计算current_row效率current_col非常低。（请记住，我们仍在谈论数十亿次评估。）
3.对P的所有元素进行操作，然后删除元素
我的另一个想法是独立计算所有有效和无效的组合。但不幸的是，由于求和错误，以下陈述是正确的：
calc_non_skipped() != calc_all() - calc_skipped()
是否有一种方便的、已知的、高性能的方法可以从初始值中获得所需的结果？
我知道这个问题相当复杂，而且相关性可能有限。尽管如此，我希望一些启发性的答案能帮助我解决我的问题。
目前的实施
目前，这通过 OpenMP 实现为 CPU 代码。我首先设置了一个上述combos 的向量，存储需要计算的每个 P 并将其传递给并行 for 循环。每个线程都有一个私有结果向量，并行区域末尾的临界区用于适当的求和。