网格三角形上的复杂OpenMP并行化

所以我有一个网格，对于每个三角形，我需要计算一些信息——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内存。