25

В качестве ответа на вопрос Какой алгоритм использовать для решения задачи? написал свою программку (приведена ниже) - ветвления с отсечением.
Бывает, меня, как поющего Кобзона :), не остановить - словом, мне захотелось выжать все, что можно. В однопоточном режиме, как мне кажется, выжал все, что мог (но, конечно, дальнейшее ускорение приветствуется). Захотел попробовать ускорить за счет параллельности - и вот тут я застрял. Любой из вариантов у меня оказывался резко хуже однопоточного.

Начал я с того, что создавал потоков по количеству ядер. Каждый поток поочередно (с синхронизацией с помощью мьютексов) брал очередную ветвь на первом уровне и полностью обрабатывал ее. Получилось плохо, как как каждая ветвь отрабатывалась полностью независимо, т.е. искалось полное решение для нее - в то время как она, быть может, была бы отброшена сразу - из-за наличия другой более короткой ветви.

Далее сделал общим значение достигнутого минимума среди всех потоков - но это значение опять же пришлось защищать мьютексом, и это привело к очередному увеличению времени работы.

Посему хотелось бы посмотреть на то, как эту задачу решат реальные эксперты в области многопоточности (к каковым себя отнести ну никак не могу).

Так сказать, конкурс на самый быстрый вариант. Интересуют именно параллельные вычисления.

А вот обещанный однопоточный код (да, я знаю, что он написан ужасно - обсуждается не это, ладно?):

#include <vector>
#include <iostream>
#include <iomanip>

using namespace std;

struct pnt { double x, y; };

using distance_t = vector<vector<double>>;

struct Func
{
    long long calls;
    double min;              // Текущее минимальное расстояние
    vector<int>        save; // Сохраненная перестановка
    const distance_t  &dist; // Расстония между точками

    Func(const distance_t& dist, const vector<pnt>& x):dist(dist)
    {
        calls = 0;
        min = 0.0;
        // Инициализация путем 0-1-2-...
        save.push_back(0);
        for(size_t i = 1; i < x.size(); ++i)
        {
            min += dist[i-1][i];
            save.push_back(i);
        }
    }

    // Проверка ветви
    bool operator()(const vector<int>& x, size_t l, double cur)
    {
        ++calls;
        // Режем все, где неверный конец (неверное начало невозможно)
        if (l != x.size()-1 && x[l] == int(x.size() - 1)) return false;

        // Текущая длина + расстояние до последней точки, если еще не достали
        if (l != x.size()-1) cur += dist[x[l]][x.size()-1];

        // Если больше минимальной - режем эту ветвь
        if (cur > min + min*DBL_EPSILON) return false;
        // Сохранение нового пути
        if (l == x.size()-1)
        {
            if (abs(min-cur) < min*DBL_EPSILON)
            {
                // Только лексикографически меньший путь
                if (save > x) save = x;
            }
            else
            {
                min = cur;
                save = x;
            }
        }
        return true;
    }
};

// Ветвление и обрезка
bool branches(size_t N, Func& f, const distance_t& dist, size_t level = 1,
              vector<int>*v_ = nullptr, double cur_dist = 0.0)
{
    // Вспомогательный вектор перестановок пути
    vector<int> * vv = (level == 1) ? new vector<int> : v_;
    if (level == 1) for(size_t i = 0; i < N; ++i) vv->push_back(i);
    vector<int>& v = *vv;

    for(size_t i = level; i < N; ++i)
    {
        // Очередная перестановка
        std::swap(v[level],v[i]);
        // Длина для нее
        double length = cur_dist + dist[v[level]][v[level-1]];
        if (f(v,level,length) && level < N-1) branches(N,f,dist,level+1,vv,length);
        // Возвращаем все, как было
        std::swap(v[i],v[level]);
    }
    if (level == 1) delete vv;
    return true;
}


int main(/*int argc, const char * argv[]*/)
{
    int N;
    cin >> N;
    vector<pnt> x;
    for(int i = 0; i < N; ++i)
    {
        // cout <<"Point " << (i+1) << ": ";
        double xx, yy;
        cin >> xx >> yy;
        x.push_back(pnt{xx,yy});
    }

    distance_t dist(N,vector<double>(N,0.0));
    for(int i = 0; i < N; ++i)
        for(int j = i+1; j < N; ++j)
            dist[i][j] = dist[j][i] = sqrt((x[i].x-x[j].x)*(x[i].x-x[j].x)+
                                           (x[i].y-x[j].y)*(x[i].y-x[j].y));

    Func f(dist,x);
    branches(N,f,dist);
    cout << "\n" << f.min << endl;
    for(auto i: f.save) cout << (i+1) << "  ";
    cout << endl;

    cout << "Calls Func(): " << setw(12) << f.calls << endl;

}
15
  • А профайлером смотрели, где основной затык? Кстати, для вычисления расстояние для декартовой плоскости есть функция hypot. Оу ... вижу комменты про его медлительность.
    – αλεχολυτ
    6 фев 2017 в 10:27
  • 4
    Дело в том, что сам алгоритм должен быть параллельным, то есть это должно быть другое решение.
    – Cerbo
    10 фев 2017 в 7:59
  • 3
    Насколько я понимаю, в мультипоточности важен принцип существования независимых контекстов - т.е имеется определенная работа которую нужно сделать, не зависящую от результатов работы в других потоках. Если это не удовлетворимо, то начинается всякая мутотень с мъютексами и кросс-потоковыми сообщениями, семафорами и др. что в итоге создаёт неудобноваримый код, который работал бы быстрее в одном потоке, не ломая людям мозги. 11 фев 2017 в 16:56
  • 1
    Тут дело такое, я конечно специалист по параллельным вычислениям, но первое, что бросилось в глаза - крайне неоптимальный алгоритм самого решения. Вот все пишут, что "NP", "перебор", все дела. На самом деле задача может быть решена гораздо быстрее методом динамического программирования по подмножествам. И работать это будет куда быстрее даже на одном ядре какого-нибудь древнего калькулятора, чем вот при таком параллельном алгоритме. Вот статья, где есть пример задачи, похожей на вашу: habr.com/ru/post/191498 Вот как раз такой алгоритм и параллелить проще. Мало ли, пригодится вам.
    – Zealint
    7 ноя 2019 в 14:38
  • 2
    По оптимизации... Мда.. dist[i][j] = dist[j][i] = ... а ты тут циклы в половиеу сократить не хочешь, если идёшь по диагонали?(сам догадаешься думаю как... обычно вообще не заполняется одна из частей матрицы симметричной) Не юзай пушбэки, если знаешь размер. Много всякого ещё, но тебя это видимо не интересует... По паралельности... Не пользуй мьютексы, пользуй атомики. А вообще Потоки, Ядра и прочее, это из пушки по воробьям. Циклы лучше паралель через openMP.
    – YepThatsMe
    9 дек 2019 в 9:42

1 ответ 1

1

Вот как надо:

#include <vector>
#include <iostream>
#include <iomanip>
#include <math.h>

using namespace std;

#define DBL_EPSILON 0.0000000000001

struct pnt { double x, y; };

using distance_t = vector<vector<double>>;

struct Func
{
    long long calls;
    double min;              // Текущее минимальное расстояние
    vector<int>        save; // Сохраненная перестановка
    const distance_t  &dist; // Расстония между точками

    Func(const distance_t& dist, const vector<pnt>& x):dist(dist)
    {
        calls = 0;
        min = 0.0;
        // Инициализация путем 0-1-2-...
        save.push_back(0);
        for(size_t i = 1; i < x.size(); ++i)
        {
            min += dist[i-1][i];
            save.push_back(i);
        }
    }

    // Проверка ветви
    bool operator()(const vector<int>& x, size_t l, double cur)
    {
        ++calls;
        // Режем все, где неверный конец (неверное начало невозможно)
        if (l != x.size()-1 && x[l] == int(x.size() - 1)) return false;

        // Текущая длина + расстояние до последней точки, если еще не достали
        if (l != x.size()-1) cur += dist[x[l]][x.size()-1];

        // Если больше минимальной - режем эту ветвь
        if (cur > min + min*DBL_EPSILON) return false;
        // Сохранение нового пути
        if (l == x.size()-1)
        {
            if (abs(min-cur) < min*DBL_EPSILON)
            {
                // Только лексикографически меньший путь
                if (save > x) save = x;
            }
            else
            {
                min = cur;
                save = x;
            }
        }
        return true;
    }
};

// Ветвление и обрезка
bool branches(size_t N, Func& f, const distance_t& dist, size_t level = 1,
              vector<int>*v_ = nullptr, double cur_dist = 0.0)
{
    // Вспомогательный вектор перестановок пути
    vector<int> * vv = (level == 1) ? new vector<int> : v_;
    if (level == 1) for(size_t i = 0; i < N; ++i) vv->push_back(i);
    vector<int>& v = *vv;

    for(size_t i = level; i < N; ++i)
    {
        // Очередная перестановка
        std::swap(v[level],v[i]);
        // Длина для нее
        double length = cur_dist + dist[v[level]][v[level-1]];
        if (f(v,level,length) && level < N-1) branches(N,f,dist,level+1,vv,length);
        // Возвращаем все, как было
        std::swap(v[i],v[level]);
    }
    if (level == 1) delete vv;
    return true;
}


int main(/*int argc, const char * argv[]*/)
{
    int N;
    cin >> N;
    vector<pnt> x;
    for(int i = 0; i < N; ++i)
    {
        // cout <<"Point " << (i+1) << ": ";
        double xx, yy;
        cin >> xx >> yy;
        x.push_back(pnt{xx,yy});
    }

    distance_t dist(N,vector<double>(N,0.0));
    for(int i = 0; i < N; ++i)
        for(int j = i+1; j < N; ++j)
            dist[i][j] = dist[j][i] = sqrt((x[i].x-x[j].x)*(x[i].x-x[j].x)+
                                           (x[i].y-x[j].y)*(x[i].y-x[j].y));

    Func f(dist,x);
    branches(N,f,dist);
    cout << "\n" << f.min << endl;
    for(auto i: f.save) cout << (i+1) << "  ";
    cout << endl;

    cout << "Calls Func(): " << setw(12) << f.calls << endl;

}

А вот ошибки: надо #include <math.h> добавить и задать DBL_EPSILON.

1
  • 1
    Ой. Я уже и забыл, это было так давно... Надо сесть, восстановить и посмотреть. В любом случае спасибо :)
    – Harry
    30 мар 2022 в 7:55

Ваш ответ

Нажимая на кнопку «Отправить ответ», вы соглашаетесь с нашими условиями использования и подтверждаете, что прочитали и поняли наши политику конфиденциальности и нормы поведения.

Всё ещё ищете ответ? Посмотрите другие вопросы с метками или задайте свой вопрос.