Условия гонки OpenMP при поиске ближайшей пары

Я выполняю задание, чтобы найти ближайшую пару между двумя непересекающимися множествами A и B. Я использую OpenMP для распараллеливания рекурсии алгоритма, но я сталкиваюсь с некоторыми гонками данных. Я очень новичок в OpenMP, поэтому я думаю, что это как-то связано с неправильным приватизацией / разделением переменных. Я поставил полный алгоритм ниже:

float OMPParticleSim::efficient_closest_pair(int n, vector<Particle> & p, vector<Particle> & q)
{
// brute force
if(n <= 3) {
float m = numeric_limits<float>::max();

for(int i = 0; i < n - 2; i++) {
for(int j = i + 1; j < n - 1; j++) {
if((set_A.find(p[i].id) != set_A.end() && set_A.find(p[j].id) != set_A.end()) || (set_B.find(p[i].id) != set_B.end() && set_B.find(p[j].id) != set_B.end())) {
continue;
}

float distsq = pow(p[i].x - p[j].x, 2) + pow(p[i].y - p[j].y, 2) + pow(p[i].z - p[j].z, 2);
pair<pair<Particle, Particle>, float> pa = make_pair(make_pair(p[i], p[j]), sqrt(distsq));

#pragma omp critical
insert(pa);

m = min(m, distsq);
}
}

return sqrt(m);
}

// copy first ceil(n/2) points of p to pl
vector<Particle> pl;
int ceiling = ceil(n/2);

for(int i = 0; i < ceiling; i++) {
pl.push_back(p[i]);
}

// copy first ceil(n/2) points of q to ql
vector<Particle> ql;

for(int i = 0; i < ceiling; i++) {
ql.push_back(q[i]);
}

// copy remaining floor(n/2) points of p to pr
vector<Particle> pr;

for(int i = ceiling; i < p.size(); i++) {
pr.push_back(p[i]);
}

// copy remaining floor(n/2) points of q to qr
vector<Particle> qr;

for(int i = ceiling; i < q.size(); i++) {
qr.push_back(p[i]);
}

float dl, dr, d;

#pragma omp task firstprivate(pl, ql, p, q, n) private(dl) shared(closest_pairs)
dl = efficient_closest_pair(ceil(n / 2), pl, ql);

#pragma omp task firstprivate(pl, ql, p, q, n) private(dr) shared(closest_pairs)
dr = efficient_closest_pair(ceil(n / 2), pr, qr);

#pragma omp taskwait
d = min(dl, dr);

float m = p[ceil(n / 2) - 1].x;
vector<Particle> s;

for(int i = 0; i < q.size(); i++) {
if(fabs(q[i].x - m) < d) {
s.push_back(Particle(q[i]));
}
}

int num = s.size();
float dminsq = d * d;

for (int i = 0; i < num - 2; i++) {
int k = i + 1;

while(k <= num - 1 && pow(s[k].y - s[i].y, 2) < dminsq) {
if((set_A.find(s[i].id) != set_A.end() && set_A.find(s[k].id) != set_A.end()) || (set_B.find(s[i].id) != set_B.end() && set_B.find(s[k].id) != set_B.end())) {
k++;
continue;
}

float dist = pow(s[k].x - s[i].x, 2) + pow(s[k].y - s[i].y, 2) + pow(s[k].z - s[i].z, 2);
pair<pair<Particle, Particle>, float> pa = make_pair(make_pair(s[i], s[k]), sqrt(dist));

#pragma omp critical
insert(pa);

dminsq = min(dist, dminsq);
k++;
}
}

return sqrt(dminsq);
}

вставить Метод выглядит так:

void OMPParticleSim::insert(pair<pair<Particle, Particle>, float> & pair) {
if(closest_pairs.size() == 0) {
closest_pairs.push_back(pair);
return;
}

for(int i = 0; i < closest_pairs.size(); ++i) {
if(closest_pairs[i].second > pair.second) {
closest_pairs.insert(closest_pairs.begin() + i, 1, pair);
break;
}
}

if(closest_pairs.size() > k) {
closest_pairs.pop_back();
}
}

Начало параллельной области здесь:

void OMPParticleSim::do_closest_pair(int num_threads) {
vector<Particle> p = set;

// presort on x
sort(p.begin(), p.end(), sortxomp);

vector<Particle> q = p;

// presort on y
sort(q.begin(), q.end(), sortyomp);
float cp;

#pragma omp parallel num_threads(num_threads)
{
#pragma omp single
{
cp = efficient_closest_pair(set.size(), p, q);
}
}

sort(closest_pairs.begin(), closest_pairs.end(), sortpairsomp);
}

Все результаты хранятся в списке closest_pairs и вывод в файл. Причина, по которой я знаю, что существуют гонки данных, заключается в том, что некоторые идентификаторы частиц имеют отрицательные значения (все они начинаются с положительных значений), и многократный запуск программы приводит к записи в файл разных значений. Любая помощь будет отличной!