一个二次算法的并行化

Question

假设我有一个包含N个元素的向量（数组，列表等），假设为V₀到V_(N-1)。对于每个元素V_i，需要计算每个索引j（包括i=j）的函数f(V_i,V_j)。该函数是对称的，因此一旦计算出f(V_i,V_j)，就无需重新计算f(V_j,V_i)。然后总共有N(N+1)/2次函数求值，使得这是一个O(N²)的算法。假设计算f所需的时间相对较长但是一致的。

现在，我想并行执行算法。我需要确定（某些数量M的）工作线程的调度，以便两个线程不同时使用内存的相同部分（即相同的元素）。例如，可以并行评估f(V₁,V₂)和f(V₃,V₄)，但不能并行评估f(V₂,V₃)。工作流程被分成步骤，每个步骤中，每个工作线程执行一次f的评估。然后对线程进行同步，然后它们继续进行下一步。

问题是，如何确定每个线程的调度（最好是最优的），作为索引对(i,j)的一系列，以便解决完整的问题（即考虑对称性访问每个索引对恰好一次）？ 尽管直接回答当然很好，但我也会感激一个指向算法或甚至相关网站/文献的指针。

- Jsl

---

为什么不能同时计算f(V1,V2)和f(V2,V3)？这个算法是否会原地修改向量？ - Alex F

并非如此，但函数f（Vi，Vj）会写入与V相同长度的另一个向量，在索引i和j处进行写入，这些写入需要受到保护。也许我应该明确说明这一点。 - Jsl

f() 函数与第二个向量做了什么？如果可以进行调整以允许同时处理重叠的成对数据，则分配任务的问题变得容易得多。例如，第二个向量是否用于记录结果？如果是这样，f() 是否可以改为将结果写入到 2D 向量的 [i][j]（和可能是[j][i]）- 然后一旦所有并行工作完成，您只需要将该 2D 向量合并到所需的 1D 版本即可。 - DMA57361

问题太笼统了，你正在要求将你的黑盒并行化。如果结果向量中的值被覆盖而不是（例如与前一个值相加），那么您只需要首先进行N次计算即可。 (v[x]，v[x+1]...v[x+1...n]) 将是无用的，因为这些计算的结果都将被随后的结果覆盖。 - Louis Ricci

Answer 1

这让我想起了一个常见的体育赛程问题：在一个有N个队的联赛中，安排N-1个比赛日，使得每个队每天都有一场比赛，并与其他所有队各打一次。

在下棋时，这个问题有一个非常形象的解决方案。将所有棋盘排成一排放在长桌上。一个玩家始终坐在同一把椅子上。其他玩家顺时针围着桌子旋转，跳过那个玩家。

- catalyst

---

太好了！我想肯定有类似的问题存在。只是补充一下（虽然我猜你已经意识到了），你应该将此方案用于i！= j的情况，然后再单独处理i = j的情况。 - Jsl

有趣！我相信这只适用于N为偶数的情况。如果N是奇数，您可以旋转“座位”，其中一个座位是休息点，没有固定的座位。 - Mikeb

Answer 2

我看到两种可能性。 一种是事先将K = N（N + 1）/ 2个任务分配给M个线程。以下只是伪代码：

allocateTasks() {
    num = ceil(N*(N+1)/2/M); // number of tasks per thread
    i0 = 0, j0 = 0;
    n = 0; 
    for (i = 0; i < N; ++i) {
        for (j = i; j < N; ++j) { // skip symmetries
            if (n == 0) {i0 = i; j0 = j;}
            ++n;
            if (n == num) {
                thread(i0, j0, num); 
                n = 0; 
            }
        }
    }
    if (n > 0) thread(i0, j0, n); 
}

thread(i0, j0, num) {
    // perform num tasks starting at (i0, j0)
    i = i0; j = j0;
    for (n = 0; n < num; ++n) {
        // f[i, j] = ...
        ++j;
        if (j >= N) {++i; j = i;}
    }
}

这里的问题是，当计算f(i,j)的时间不同时，负载在线程之间就不会平衡。

因此，也许第二种方法更好。每个线程都会执行下一个可用的任务，而没有全局分配任务的预先设定。

// init boolean array: done[0..N-1, 0..N-1] = false
thread() {
    for (i = 0; i < N; ++i) {
        for (j = i; j < N; ++j) { // skip symmetries
            boolean mytask = false;
            lock(mutex); // synchronize access to global array done[][]
            if (!done[i, j]) {
                done[i, j] = true;
                mytask = true;
            }
            unlock(mutex);
            if (mytask) {
                f(i, j) = ...
            }
        }
     }
}

无论如何，访问f(i, j)的方式都是通过。

g(i, j) = j >= i ? f(i, j) : f(j, i)

- Jiri Kriz

Answer 3

让我们看一下直接实现：

for(i=0; i < N; ++i)
{
    for(j=1; j < i; ++J)
    {
        计算 i,j 对并分配给 i,j 和 j,i 结果
    }
}

我是C++程序员，所以我可以考虑使用OpenMP进行外部循环并行化：

#pragma OMP parallel for
for(i=0; i < N; ++i)
{
    ....
}

如果您不了解OpenMP，它只是将循环分成n个循环，根据处理器的数量，并在并行执行它们。我认为在这种情况下结果不会很好，因为每个i+1循环都比i循环短。让我们按照这种方式编写此算法，使所有循环具有相同的长度。总循环次数将为N/2。第一个循环处理0和N-1行。第二个循环处理1和N-2行等。这样的循环可以成功并行化，而不会发生冲突。简化代码，没有关于偶数/奇数N等的细节：

#pragma OMP parallel for
for(i=0; i < N/2; ++i)
{
    处理 i 值
    ...
处理 N - 1 - i 值
    ...
}

这只是一个大致的想法，你可以修正其中的错误细节。

- Alex F