3

我正在并行化我的粒子单元代码，我用它来模拟地球内的 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];
}

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

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

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

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

有没有人有过这个或类似问题的经验？

4

2 回答 2

2

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

#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循环内的条件决定哪个线程将更新数组pos的位置。grid这个位置只属于一个线程，因此atomic不需要保护。
  3. 该公式(pos/factor)%nthreads只是一个简单的启发式。任何pos返回该范围内的值的函数0,...,nthreads-1实际上都可以替换为该表达式，而不会影响最终结果的有效性（因此，如果您有更好的选择，请随时更改它）。请注意，此功能选择不当可能会导致负载平衡问题。
于 2012-11-07T21:24:07.723 回答
1

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

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

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

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

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

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

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

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

您可以尝试来自 openMP的动态 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
      }
   }
}
于 2012-11-07T23:42:19.993 回答