爲什麼不的OpenMP並行VTK IntersectWithLine代碼

Question

我很難理解爲什麼for循環，使用OpenMP並行化，不若VTK代碼在循環內執行使用所有n_threads線程（= 2X #cores）。具體而言，我想用網格相交線條/光線。我通過

1. 建設OBB樹從網格
2. 所有需要的線相交，其中網

因爲我想並行它遵循this tutorial，我創建n_threads樹，使得每個線程可以使用它自己的樹的實例：

// Pre-allocate the array 
int n_threads = omp_get_max_threads(); 
trees = std::vector<vtkSmartPointer<vtkOBBTree>>((unsigned int) n_threads); 

// Build n_threads OBB trees 
#pragma omp parallel for num_threads(n_threads) 
for (int t = 0; t < n_threads; ++t) 
{ 
    trees[t] = vtkSmartPointer<vtkOBBTree>::New(); 
    vtkSmartPointer<vtkPolyData> mesh_clone = vtkSmartPointer<vtkPolyData>::New(); 
    #pragma omp critical (build_mesh_tree) 
    { 
     mesh_clone->DeepCopy(mesh); 
    } 
    trees[t]->SetDataSet(mesh_clone); 
    trees[t]->BuildLocator(); 
} 

然後我遍歷所有點來計算origin和在points

#pragma omp parallel for num_threads(n_threads) 
for (unsigned long i = 0; i < n_points; ++i) 
{ 
    vtkSmartPointer<vtkPoints> intersection_points = vtkSmartPointer<vtkPoints>::New(); 
    vtkSmartPointer<vtkIdList> cell_ids = vtkSmartPointer<vtkIdList>::New(); 

    int this_thread = omp_get_thread_num(); 
    int code = trees[this_thread]->IntersectWithLine(
      origin.data(),  // pointer to the raw data 
      points.at(i).data(), // buffer of a Eigen matrix 
      intersection_points, 
      cell_ids); 

    // Do something with intersection_points and cell_ids 
} 

OpenMP的各點表明如預期的簡單的C++代碼工作。但是，當纏繞Vtk呼叫時，它不能達到目的。我想這是由於Vtk已經提供了一個parallelization framework（ref. to the guide）。

如果是這樣，你能解釋一下，爲什麼OpenMP無法並行運行vtk相關的代碼？如果不是，可能是什麼原因？

Answer 1

它究竟如何失敗？你剛剛嘗試過嗎？打印線程號以查看是否生成了n_threads？如果您只是在intersection_points和cell_ids中得到錯誤結果，那是因爲每當它們進入IntersectWithLine函數（您可以查看執行here，第800-808行）時，每個線程都會重置這些陣列。

爲了解決這個問題，想到的最簡單的解決方案是每個線程都有自己的列表，然後在關鍵部分將它們連接起來，但是更快可能會爲每個線程預先分配這些列表的數組，然後閱讀從它的結果，例如：

vtkSmartPointer<vtkPoints> *inter_points = new vtkSmartPointer<vtkPoints> [n_threads]; 
vtkSmartPointer<vtkIdList> *inter_cells = new vtkSmartPointer<vtkIdList> [n_threads]; 
for (unsigned long i = 0; i < n_threads; ++i) 
{ 
    inter_points[i] = vtkSmartPointer<vtkPoints>::New(); 
    inter_cells[i] = vtkSmartPointer<vtkIdList>::New(); 
} 

#pragma omp parallel for num_threads(n_threads) 
for (unsigned long i = 0; i < n_points; ++i) 
{ 
    int this_thread = omp_get_thread_num(); 
    int code = trees[this_thread]->IntersectWithLine(
      origin.data(), // pointer to the raw data 
      points.at(i).data(), // buffer of a Eigen matrix 
      inter_points[this_thread], 
      inter_cells[this_thread); 
} 

// now concatenate the results to one array if necessarry 

（沒有編制，可能有語法錯誤，只是普通的C，所以不是最安全的方式......）