Reputation: 55
The problem that I am trying to solve is the following:
I have a big point set which I try to split into smaller ones. I know the sizes of the smaller ones, and I am trying to populate the smaller ones using the bigger one in parallel. In order to do that, I have an array of indices per smaller pointset, which will help me populate my point sets correctly. I have developed the following code. But the omp atomic capture
directive does not help me get the same results as I do sequentially.
// points
std::vector<Eigen::Matrix3Xd> points(numberofPointSets);
// initialize arrays
for (auto& pointSetInformation: pointSetsInformation)
{
int pointSetSize = getPointSetSize();
int index = getIndexOfPointsSet();
points[index].resize(3, pointSetSize);
}
// currentIndices to assign specific columns of the matrix in parallel
std::vector<int> currentIndices(numberofPointSets, -1);
// Assign arrays
#pragma omp parallel num_threads(numberOfThreads) shared(pointsArray, points, currentIndices) default(none)
{
int index;
int currentIndexInCurrentPointSet;
double point[3];
#pragma omp for schedule(static)
for (int i = 0; i < input->GetNumberOfPoints(); ++i)
{
index = getIndexOfPointsSet();
// update currentIndices (we start from -1, that's why this happens first)
#pragma omp atomic capture
currentIndexInCurrentPointSet = ++currentIndices[index];
pointsArray->GetPoint(i, point);
points[index](0, currentIndexInCurrentPointSet) = point[0];
points[index](1, currentIndexInCurrentPointSet) = point[1];
points[index](2, currentIndexInCurrentPointSet) = point[2];
}
}
Am I missing something? Can the omp atomic capture
directive update array cells, or only variables?
Thanks a lot in advance.
EDIT:
Is there a way to enforce the same order of operation as the i
variable?
Upvotes: 1
Views: 535
Reputation: 51393
Is there a way to enforce the same order of operation as the i variable?
You try to use the OpenMP ordered clause
#pragma omp parallel num_threads(numberOfThreads) shared(pointsArray, points, currentIndices) default(none)
{
int index;
int currentIndexInCurrentPointSet;
double point[3];
#pragma omp for ordered schedule(static, 1)
for (int i = 0; i < input->GetNumberOfPoints(); ++i)
{
#pragma omp ordered
{
index = getIndexOfPointsSet();
currentIndexInCurrentPointSet = ++currentIndices[index];
}
pointsArray->GetPoint(i, point);
points[index](0, currentIndexInCurrentPointSet) = point[0];
points[index](1, currentIndexInCurrentPointSet) = point[1];
points[index](2, currentIndexInCurrentPointSet) = point[2];
}
}
Upvotes: 2