从随机分布的粒子到规则网格的通信优化并行化

Question

我正在并行化我的粒子单元代码，我用它来模拟地球内部的 2D 和 3D 变形。使用 OpenMP 可以轻松并行化代码的几个例程，并且可以很好地扩展。但是，我在处理从粒子到网格单元的插值的代码的关键部分中遇到了问题。粒子在每次迭代中四处移动（根据速度场）。许多计算在规则的、不变形的网格上执行是最有效的。因此，每次迭代都涉及从“随机”分布的粒子到网格单元的通信。

这个问题可以用以下简化的一维代码来说明：

//EXPLANATION OF VARIABLES (all previously allocated and initialized, 1D arrays)
//double *markerval; // Size Nm. Particle values. Are to be interpolated to the grid
//double *grid; // Size Ng=Nm/100 Grid values. 
//uint *markerpos; // Size Nm. Position of particles relative to grid (each particle
// knows what grid cell it belongs to) possible values are 0,1,...Ng-1

//#pragma omp parallel for schedule(static) private(e)
for (e=0; e<Nm; e++) {
    //#pragma omp atomic
    grid[markerpos[e]]+=markerval[e];
}

在最坏的情况下粒子位置是随机的，但通常情况下，粒子在内存中彼此相邻，在空间中也彼此相邻，因此也在网格内存中。

如何有效地并行化这个过程？几个粒子映射到同一个网格单元，因此如果上述循环直接并行化，竞争条件和缓存交换的可能性非零。使更新原子化可以防止竞争条件，但会使代码比顺序情况慢得多。

我还尝试为每个线程制作一个网格值的私有副本，然后随后将它们相加。然而，这可能需要在代码中使用太多内存，并且对于这个示例，它并不能很好地随线程数量而扩展（原因我不确定）。

第三种选择可能是从网格映射到粒子，然后循环通过网格索引而不是粒子索引。然而，我担心这会涉及很多，并且需要对代码进行重大更改，而且我不确定它会有多大帮助，因为它还需要使用计算成本也很高的排序例程。

有没有人遇到过这个或类似的问题？

Answer 1

一个选项可以是在线程上手动映射迭代：

#pragma omp parallel shared(Nm,Ng,markerval,markerpos,grid)
{
  int nthreads = omp_get_num_threads();
  int rank     = omp_get_thread_num();
  int factor   = Ng/nthreads;

  for (int e = 0; e < Nm; e++) {
    int pos = markerpos[e];
    if ( (pos/factor)%nthreads == rank )
      grid[pos]+=markerval[e];
  }
}

几点说明：

1. for 循环的迭代不在线程之间共享。而是每个线程执行所有迭代。
2. for 循环内的条件决定哪个线程 将更新grid 数组的位置pos。此位置只属于一个线程，因此不需要atomic 保护。
3. 公式(pos/factor)%nthreads 只是一个简单的启发式。 pos 的任何返回值在0,...,nthreads-1 范围内的函数实际上都可以替换为该表达式，而不会影响最终结果的有效性（因此，如果您有更好的机会，请随时更改它）。请注意，此功能选择不当可能会导致负载平衡问题。

Answer 2

我还使用 OpenMP 并行化了分子动态算法。首先，您必须分析算法瓶颈（例如内存限制和 CPU 限制）。这样你就知道哪里需要改进了。

最初，我的 MD 受内存限制，因此只需将数据布局从结构数组 (AOS) 更改为数组结构 (SOA)，我的速度大约提高了 2x（由于空间局部性）。对于只适合 RAM 的输入，我还应用了一种阻塞技术。原始算法计算每个粒子之间的力对如下：

for(int particleI = 0; i < SIZE ; i++)
 for(int particleJ = 0; j < SIZE; j++)
     calculate_force_between(i,j);

基本上，使用块技术，我们通过粒子块来聚合力计算。例如，计算前 10 个粒子之间的所有力比，然后计算下 10 个粒子，以此类推。

这种blocks技术的使用促进了对时间局部性的更好利用，因为使用这种方法可以在更短的时间内实现对相同粒子的更多计算时间。因此，降低了我们尝试访问的值不再在缓存中的可能性。

现在我有一个 MD CPU 绑定，我可以尝试使用 multi-threads 来改进它，但首先，你需要：

1. 验证您的算法在哪里花费了大部分执行时间；
2. 找到可以并行完成的任务并确定它们的 粒度（检查其并行化是否合理）；
3. 负载均衡，确保线程间工作负载均衡；
4. 尽量减少同步的使用。

由于负载平衡问题，我在扩展我的 MD 时遇到了问题。一些线程比其他线程做更多的工作。解决方案？

您可以从 openMP 尝试 dynamic for。请注意，在 OpenMP 中，您可以指定要分配给线程的工作块。但是，在定义块时必须小心！使用动态 for，块太小会导致同步开销，太大会导致负载平衡问题。

我也遇到了同步开销问题。我使用的是关键算法，但算法无法扩展。我用更细粒度的同步替换了这个关键，即锁，每个粒子一个。我对这种方法进行了一些改进。

作为最后一种方法（处理同步开销），我使用数据冗余。每个粒子完成它的工作并将结果保存在一个私有的临时数据结构中。最后，所有线程都降低了它们的值。在所有版本中，这是给我最好结果的版本。

我能够在 CPU 上实现良好的加速，但与我在 GPU 版本中实现的速度相比却没有。

根据您提供的信息，我会这样做：

omp_lock_t locks [grid_size]; // create an array of locks
int g;
#pragma omp parallel for schedule(static)
for (e=0; e<Nm; e++)
{
    g = markerpos[e];

    omp_set_lock(&locks[g]);
    grid[g]+=markerval[e];
    omp_unset_lock(&locks[g]);
}

从，我理解的问题是你必须使用atomic 来确保多个线程不会同时访问同一个抓握位置。作为一种可能的解决方案，您可以创建一组锁，并且每次线程必须访问网格的一个位置时，它都会请求并获取与该位置关联的锁。另一种解决方案可以是：

double grid_thread[grid_size][N_threads]; // each thread have a grid
// initialize the grid_threads to zeros

#pragma omp parallel
{
    int idT = omp_get_thread_num();
    int sum;
    #pragma omp parallel for schedule(static)
    for (e=0; e<Nm; e++)
        grid_thread[markerpos[e]][idT]+=markerval[e]; // each thread compute in their 
                                                     // position
    for(int j = 0; j <Nm; j++)
    { 
        sum = 0;
        #pragma omp for reduction(+:sum) 
        for (i = 0; i < idT; i++)                   // Store the result from all
           sum += grid_thread[j][i];                // threads for grid position j

         #pragma barrier                            // Ensure mutual exclusion

         #pragma master
         grid[j] +=sum;                             // thread master save the result  
                                                    // original grid
         #pragma barrier                            // Ensure mutual exclusion
      }
   }
}