我正在处理与我的粒子中的代码并行，我用来在2D和3D中进行地球内变形的模拟。代码的几个例程使用OpenMP易于并行化，并且可以很好地缩放。但是，我在代码的关键部分中遇到了问题，该代码涉及从粒子到网格单元的插值。颗粒正在为每次迭代（根据速度场）移动。许多计算最有效地在常规的，不形成的网格上执行。因此，每次迭代都涉及从"随机"分布的粒子到网格单元的通信。
可以在以下简化的1D代码中说明问题：
//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];
}
在最坏的情况下，粒子位置是随机的，但通常是在记忆中相邻的粒子，也在空间中相邻邻居，因此在网格记忆中也相邻。
如何有效地平行此过程？几个粒子映射到同一网格单元格，因此，如果上述环直接并行化，则种族条件的机会非零的机会交换。使更新原子可以防止种族条件，但使代码比顺序案例慢得多。
我还试图为每个线程制作一个私人副本，然后随后添加它们。但是，这可能需要在代码中使用过多的内存，并且在此示例中，线程数量的扩展不是很好（出于我不确定的原因）。
第三个选项可能是从网格映射到粒子，然后通过网格指数而不是粒子指数循环。但是，我担心，这将非常参与其中，并且需要对代码进行重大更改，而且我不确定它会有所帮助，因为它也需要使用排序例程，而排序例程也将在计算上也很昂贵。p>有人在此或类似问题上有任何经验吗？
一个选项可以是手动在线程上映射迭代：
#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];
}
}
一些言论：
-
for循环的迭代在线程之间未共享。而是每个线程进行所有迭代。 -
for循环内部的条件性决定哪个线程将更新grid数组的位置pos。此位置仅属于一个线程，因此不需要atomic保护。 - 公式
(pos/factor)%nthreads只是一个简单的启发式。pos的任何函数返回在0,...,nthreads-1范围内的值的任何函数实际上可以代替此表达式，而不会损害最终结果的有效性（因此，如果您有更好的镜头，请随时更改它）。请注意，此功能的不良选择可能导致负载平衡问题。
我也与OpenMP平行了分子动力学算法。首先，您必须分析算法瓶颈（，例如，内存绑定和CPU绑定）。这样，您将知道在哪里改进。
最初，我的MD是内存绑定的，因此我仅通过将数据布局从结构（AOS）更改为阵列（SOA）（SOA）（适当的空间locatial ）。我还申请了仅适用于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来改进它，但是首先，您需要：
- 验证您的算法在何处花费其大部分执行时间;
- 找到可以并行完成的任务并确定其粒度（检查其并行化是否合理）;
- 加载余额，确保线程之间的工作良好载荷;
- 最小化同步的使用。
由于负载平衡问题，我在扩展MD方面遇到了问题。有些线程比其他线程更多的工作。解决方案？
您可以从OpenMP尝试的动态。请注意，在OpenMP中，您可以指定要分配给线程的工作块。但是，您必须小心定义块！对于的动态，块太小会导致开销的同步，太大会导致负载平衡问题。
我也有同步开销的问题。我正在使用关键，算法没有扩展。我用较细的谷物同步代替了关键，即锁，每个粒子一个。我对这种方法有所改进。
作为最后一种方法（处理同步开销），我使用数据冗余。每个粒子都进行了工作，并将结果保存在私人临时数据结构中。最后，所有线程都降低了其值。从所有版本中，这是给我最好的结果的版本。
我能够在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
}
}
}