動的計画法を使用して最小コストのポリゴン三角形分割を計算し、OpenMPを使用してそれを並列化/ベクトル化する関数を作成しようとしています。これまでに作成したコードは正しい結果を返しますが、速度が遅すぎます。3000ポイントを超えるポリゴンの場合、5分経っても停止しません。コードは次のとおりです。
#pragma omp declare simd
float dist(float x1, float y1, float x2, float y2)
{
return sqrt((x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2));
}
float triangulate(const vector<Point> &points) {
int n = points.size();
vector<vector<float>> table (n, vector<float>(n, 0));
int threads = omp_get_max_threads();
for (int gap = 0; gap < n; ++gap)
{
for (int i = 0, j = gap; j < n; ++i, ++j)
{
if (j < i+2)
table[i][j] = 0.0;
else
{
int size = j - i - 1;
Point p1 = points[i], p2 = points[j];
//Precompute distance between i and j
float ij = dist(p1.x, p1.y, p2.x, p2.y);
float minimum = MAX;
#pragma omp parallel for simd schedule(static, 64) num_threads(threads) reduction(min:minimum) if(size > 300)
for (int k = i+1; k < j; ++k)
{
Point p3 = points[k];
float perimeter = ij + dist(p1.x, p1.y, p3.x, p3.y) + dist(p2.x, p2.y, p3.x, p3.y) + table[i][k] + table[k][j];
if(perimeter < minimum)
{
minimum = perimeter;
}
}
table[i][j] = minimum;
}
}
}
return table[0][n-1];
}
ギャップとi、j forループimhoは並列化できないため、k上のforループのみを並列化できます。私はスケジュールの議論と一緒に遊んでみましたが、改善はありませんでした。私は何かが足りないのですか、それともこのアプローチではこの機能だけを速くすることはできませんか？
i/で並列化してみませんjか？与えられたI / jのペアについて、境界を計算するだけの値に依存しているtable[i][k]とtable[k][j]の間のギャップiとkの間、またはkとjの間よりも小さくなるiとj。最も外側のループgapが順番に実行されている限り、i/の内側のループにjは恥ずかしい並列処理プロパティがあり、予防策なしで並列化できます。
// Do this here so that we can dispose of the if/then/else in the loop later
table[0][0] = 0.0;
for (int i= 1; i < n; ++i){
table[i][i-1] = 0.0;
table[i][i] = 0.0;
}
// spawn parallel threads here
#pragma omp parallel default(shared)
for (int gap = 2; gap < n; ++gap)
{
// Loop on i will now be distributed among threads
// Not sure that is possible to place two variables in the loop definition
// so loop on i and compute the corresponding value of j
#pragma omp for
for (int i = 0; i < n-gap; ++i)
{
int j=i+gap;
int size = gap - 1;
Point p1 = points[i], p2 = points[j];
//Precompute distance between i and j
float ij = dist(p1.x, p1.y, p2.x, p2.y);
float minimum = MAX;
// remove of all directives but simd
#pragma omp simd reduction(min:minimum)
for (int k = i+1; k < j; ++k)
{
Point p3 = points[k];
float perimeter = ij + dist(p1.x, p1.y, p3.x, p3.y) + dist(p2.x, p2.y, p3.x, p3.y) + table[i][k] + table[k][j];
if(perimeter < minimum)
minimum = perimeter;
}
table[i][j] = minimum;
} // implicit barrier, all threads will wait here for loop completion before moving to next value of gap
} //end of loop, end of parallel region
この記事はインターネットから収集されたものであり、転載の際にはソースを示してください。
侵害の場合は、連絡してください[email protected]
コメントを追加