Динамическое программирование распараллеливания / векторизации с OpenMP

Question

Я пытаюсь написать функцию, которая вычисляет минимальную стоимость триангуляции многоугольника, используя динамическое программирование, и парализует / векторизует ее, используя 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 для циклов imho нельзя распараллелить, поэтому можно распараллелить только цикл for над k.Я пытался подыгрывать графику, но без улучшений.Я что-то упустил или просто эта функция не может быть быстрее в рамках этого подхода?

Answer 1

Почему бы не распараллелить на 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