8

N 体问题的一个非常常见的问题是使用双循环来计算粒子之间的相互作用。考虑具有 n 个粒子的 N 体问题，循环可以写成

for (i = 0, i < n; i++)
    for (j = i+1, j < n; j++)
        // calculate interaction

我的问题是如何使用不同的线程并行化这个循环。目标是每个线程“理想地”必须计算相同数量的交互。

我的想法是在不同的间隔上分离外循环，即 i 循环，比如 a_k=a(k)，其中 k = 1,2,...,p 其中 p 是我们想要划分的线程数问题成。

所以，循环可以写成

for (k = 1, k < p; k++)
    for (i = a(k), i < a(k+1); i++)
        for (j = i+1, j < n; j++)
            // calculate interaction

其中最外层的循环，即 k 循环，是要并行化的循环。

因为最内层循环j-cycle的交互次数为n-(i+1)，所以每个线程计算的交互次数为

\sum_{i=a(k)}^{a(k+1)} n - (i+1)

这意味着人们想找到离散函数 a_k 使得函数

f[a_k] = \sum_{i=a(k)}^{a(k+1)} n - (i+1)

边界条件 a(1)=0 和 a(p)=n 是一个常数函数，因此迫使每个线程上的交互次数相同。

我尝试过使用不同的“启发式”（例如 a_k 多项式、指数、对数），但到目前为止还没有一个给我一个满意的答案。这个问题的直接解决方案对我来说并不明显。

对于小p，这个问题可以放在“最小化麻袋问题”上，基本上每个a_k都是一个变量来最小化函数

f(a_1,a_2,a_3,...) = 总和(|f[a_k] - n/p|^2)

但是您可能会猜到，对于较高的 p 值，这不是有效的（甚至收敛）。

有谁知道如何解决这个问题？

4

4 回答 4

3

（对不起，如果没有清楚地表达出来，这在我的脑海中是有道理的）。

将 1 到 N 的所有数字相加时，您会注意到 N + 1 = (N - 1) + 2 = (N - 2) + 3 等。

那么，如果每个线程使用一个小 i 和一个大 i，这样总和总是相加呢？

或者，假设您想始终使用 5 个线程。线程 1 将执行前 10% 和最后 10%，线程 2 将执行第二个 10% 和倒数第二个 10%，依此类推。“早期”和“晚期”部分的每一对都会加起来相同的交互总数。

编辑：

从另一个帖子中窃取图表...

   0 1 2 3 4 5 6 7 8

0  - A B C D D C B A
1    - B C D D C B A  
2      - C D D C B A
3        - D D C B A  
4          - D C B A
5            - C B A
6              - B A
7                - A
8                  -

这是否更清楚地表明了我的意思？

于 2012-05-09T23:05:04.483 回答
3

您可以将您的对象分成k几组大致的N/k物体，并使用它来将您的初始交互三角形分解为多个k*(k + 1)/2部分：

   0 1 2 3 4 5 6 7 8
                      -- N=9;  k=3;  N/k=3
0  - A A B B B C C C
1    - A B B B C C C  -- diagonal pieces:  A, D, F
2      - B B B C C C
3        - D D E E E  -- non-diagonal pieces: B, C, E
4          - D E E E
5            - E E E
6              - F F
7                - F
8                  -

由于存在两种类型的碎片，这种观点变得复杂：沿对角线的那些（带有(N/k)*(N/k - 1)/2元素的三角形）和不沿对角线的那些（带有元素的正方形(N/k)*(N/k)）。但是，由于对角线块的大小大约是方形块的一半，因此您可以为每个线程分配两个以平衡负载——总共k*k/2大致相等的任务。

这种方法的一个优点是每个任务只需要访问2*N/kbody 的数据，这可以使它对缓存更加友好。

于 2012-05-10T00:30:26.237 回答
2

假设你的编译器支持 OpenMP，你为什么不能简单地尝试做

#pragma omp parallel for schedule(dynamic) // or: schedule(guided)
for (i = 0; i < n; i++)
    for (j = i+1; j < n; j++)
        // calculate interaction

甚至（您需要进行基准测试以了解哪个表现更好）

#pragma omp parallel
const int stride = omp_get_num_threads() + 1; 
for (i = omp_get_thread_num(); i < n; i += stride)
    for (j = i+1; j < n; j++)
        // calculate interaction
于 2012-05-10T09:40:39.640 回答
0

今天我刚刚找到了解决方案。在有人确认之前我不会接受

为了使 f[a_k] 是关于 k 的常数函数，则

f[a_{k+1}] - f[a_k] = 0

k = 1,2,3,...,p-1 时必须为真。

我们可以使用我在问题上发布的定义扩展这个方程，我们得到一个关于 a_k, k=1,2,3,...,p 的“p”2º 阶代数方程系统。我没有看到任意 p 的封闭解决方案，但可以对每个 p 进行解析求解。

我已经确认：

  1. 总和，当使用我计算的 a_k 是 n(n-1)/2，这个问题的交互总数。

  2. 对于 p = 2,3,4,5 和 10，每个线程的交互次数确实是恒定的（其中 p=10 在mathematica® 上需要一些时间来计算）。

编辑

在详细检查了不同 p 值的解决方案之后，我已经达到了一般封闭解决方案

a_k = 1/(2 p) (-p + 2 pn - sqrt[p^2 + 4 p (p + 1 - k) (n - 1) n])

这对每个 p>=2, n>1 都有效。

这完成了答案。

于 2012-05-10T08:35:35.337 回答