所以我有一个网格,对于每个三角形,我需要计算一些信息——computeInfo(triangle)
。在我的顺序版本中,我只是使用迭代器来遍历整个网格。
为了获得一些时间,我正在尝试用OpenMP并行化这一部分。诀窍是computeInfo(T)
访问并修改三角形T
周围的一些三角形。
我的想法如下:首先计算网格的边界框。然后创建一个足够的三维网格(单元格足够大,这样computeInfo(T)
就不会引起任何问题),并将每个三角形分配给相应的网格单元格。
然后,每个线程处理一个单元格中的所有三角形,该单元格的坐标可以表示为(2*i, 2*j, 2*k)
:这确保了不会因另一个线程而出现修改。我们等待处理所有相应的单元格,然后处理(2*i+1, 2*j, 2*k)
、(2*i, 2*j+1, 2*k)
,依此类推,直到处理(2*i+1, 2*j+1, 2*k+1)
。
显然,我的并行代码几乎没有顺序代码快(我做的最好的事情是快2倍,使用8个线程…)。我认为这是由于我用于存储网格的结构,要处理的方面等等。
下面是我代码的简化版本:
vector<set<Index> > vsGrid(nX*nY*nZ); // vector that contains each "cells"
// a cell is a set of the facet index
// it contains
// nX, nY, nZ = size of the grid
vector<vector<set<Index> > > vvsGrid(iNbOfThreads);
// one vector for each thread
#pragma omp parallel for
for(int i = 0; i < iNbOfThreads; ++i)
{
vector<set<Index> > vsGrid(nX*nY*nZ);
vvsGrid[i] = vsGrid;
}
#pragma omp parallel
{
const int iThrdId = omp_get_thread_num();
// each thread process a part of the total triangles (called facet)
for(PlaneFinderAPI::Polyhedron::Facet_iterator f = vFStart[iThrdId]; f != vFEnd[iThrdId]; ++f)
{
Point_3d p = barycenterOf(f);
int X = int((p.x() - minX)/(12*meanEdgeSize));
int Y = int((p.y() - minY)/(12*meanEdgeSize));
int Z = int((p.z() - minZ)/(12*meanEdgeSize));
vvsGrid[iThrdId][X+nX*Y+nX*nY*Z].insert(f->index());
}
#pragma omp barrier
for(int col = iThrdId; col < nX*nY*nZ; col+=iNbOfThreads)
{
for(int j = 0; j < iNbOfThreads; ++j)
{
// we merge the cells of all the threads
vsGrid[col].insert(vvsGrid[j][col].begin(), vvsGrid[j][col].end());
}
}
}
for(int i = 0; i < vvCellToProcess.size(); ++i)
{
#pragma omp parallel for
for(int j = 0; j < vvCellToProcess[i].size(); ++j)
// vvCellToProcess contains the cells to process, first vector contains all
// the cells (2*i, 2*j, 2*k) and so on
{
const int iCellToProcess = vvCellToProcess[i][j];
for(set<Index>::iterator it = vsGrid[iCellToProcess].begin(); it != vsGrid[iCellToProcess].end(); ++it)
{
PlaneFinderAPI::Polyhedron::Facet_iterator f = facetMap.at(*it);
computeInfo(f);
}
}
}
我认为主要的问题在于过度使用集合的向量中的向量。我认为内存没有得到有效分配,这就是为什么它很慢的原因。那么我应该使用什么结构呢?有没有更有效的方法来解决我的问题?
序列号:
for (PlaneFinderAPI::Polyhedron::Facet_iterator f = P.facets_begin() ; f != P.facets_end() ; ++f)
{
computeInfo(f);
}
computeInfo
使用网格的标准半边数据结构来查看f周围的顶点。为了避免对同一个顶点进行多次查看,每个顶点都与一个布尔visited
相关联。这就是为什么我不能简单地并行化for循环,因为在读取/写入visited
布尔值时会出现问题。
技术信息:我使用的是Windows 7和Visual Studio 2010,i7 960@3.20GHz,有8个线程和12GB内存。
我建议您先对代码进行基准测试,然后再猜测什么是慢的,什么是快的。同样取决于索引是什么,每次分配和读取向量时可能会有一个昂贵的副本。您可以尝试存储索引的指针。
例如,这一行:
vvsGrid[i] = vsGrid;
将导致向量中每个集合的多个索引的副本,然后将这些集合中的每个集合和每个索引的副本复制到新向量中。