Участник:Bvdmitri/Адаптивное интегрирование

Основные авторы описания: Д.В.Багаев

1 Свойства и структура алгоритма

1.1 Общее описание алгоритма

Адаптивное интегрирование — метод быстрого вычисления интеграла функции, который использует адаптивное изменение плотности сетки интегрирования. На промежутках, где функция гладкая, сетка разреженная, на промежутках со значительными колебаниями значений сетка гуще для получения необходимой точности.

1.2 Математическое описание алгоритма

Формулы метода:

1. Интеграл функции вычисляется по квадратурной формуле Симпсона с числом узловых точек равным 4 и 8 на отрезке [math][a, b] [/math].
2. Вычисляется разность найденных значений
   1. Если разность интегралов меньше либо равна желаемой точности интегрирования ε - заканчиваем процедуру.
   2. Если разность больше, то отрезок разбивается на два : [math][a, \tfrac{a + b} {2}][/math] и [math][\tfrac{a + b} {2}, b][/math] и для каждого отрезка повторяется процедура 1, затем возвращается сумма полученных значений.

Алгоритм имеет рекурсивную природу.

1.3 Вычислительное ядро алгоритма

Основное время работы алгоритма приходится на рекурсивные вызовы, а также на вычисление определенных интегралов методом Симпсона на отрезках, полученных в результате рекурсивного разбиения.

1.5 Схема реализации последовательного алгоритма

Простейшую рекурсивную реализацию данного алгоритма , которая не использует ресурсы параллелизма можно посмотреть здесь.

#include "SequentialMethod.h"

SequentialMethod::SequentialMethod(double eps, const std::function<double(double)> &f) : IntegrationMethod(eps, f) {}

double SequentialMethod::integrate(double a, double b) {
    double x = simpson(a, b, 4);
    double y = simpson(a, b, 8);
    if (fabs(x - y) < eps) {
        return y;
    } else {
        return integrate(a, (a + b) / 2.0) + integrate((a + b) / 2.0, b);
    }
}

SequentialMethod::~SequentialMethod() {}

1.6 Последовательная сложность алгоритма

Явную формулу сложности последовательного алгоритма выписать не представляется возможным, так как сложность зависит от структуры интегрируемой функции, длины отрезка, а так же необходимой точности.

1.7 Информационный граф

У алгоритма есть несколько возможностей для распараллеливания.
В первом случае можно разбить отрезок на фиксированное число под-отрезков, и каждый под отрезок считать на одном процессоре с дальнейшим суммированием результатов с каждого процессора. Внутри процессора используется схема стека задач. Рекурсия заменяется на операции добавить и достать задачу из общего стека задач.

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

1.9 Входные и выходные данные алгоритма

1.9.1 Входные данные

• Функция.
• Отрезок интегрирования [math][a, b][/math].
• Желаемая точность интегрирования [math]\varepsilon[/math].

1.9.2 Выходные данные

• Определенный интеграл функции на отрезке [a, b] с точностью ε.

3 Литература

3.1 Исходный код реализованного алгоритма

Весь исходный код в статье можно найти в репозитории на GitHub.

3.2 Реализации на одном процессоре

Простейшую рекурсивную реализацию данного алгоритма , которая не использует ресурсы параллелизма можно посмотреть здесь.

#include "SequentialMethod.h"

SequentialMethod::SequentialMethod(double eps, const std::function<double(double)> &f) : IntegrationMethod(eps, f) {}

double SequentialMethod::integrate(double a, double b) {
    double x = simpson(a, b, 4);
    double y = simpson(a, b, 8);
    if (fabs(x - y) < eps) {
        return y;
    } else {
        return integrate(a, (a + b) / 2.0) + integrate((a + b) / 2.0, b);
    }
}

SequentialMethod::~SequentialMethod() {}

Реализация с помощью библиотеки Intel TBB для более эффективного использования многоядерных процессоров можно посмотреть здесь.

#include "TBBMethod.h"

TBBMethod::TBBMethod(double eps, const std::function<double(double)> &f) : IntegrationMethod(eps, f) {}

double TBBMethod::integrate(double a, double b) {
    double x = simpson(a, b, 4);
    double y = simpson(a, b, 8);
    if (fabs(x - y) < eps) {
        return y;
    } else {
        tbb::task_group group;

        group.run([&] {
            x = integrate(a, (a + b) / 2.0);
        });

        group.run_and_wait([&] {
            y = integrate((a + b) / 2.0, b);
        });

        return x + y;
    }
}

TBBMethod::~TBBMethod() {}

3.3 Параллельная реализации

3.3.1 Реализация с помощью MPI для межпроцессорного взаимодействия и Intel TBB для внутрипроцессорного взаимодействия

#include "MPI_TBBMethod.h"

MPI_TBBMethod::MPI_TBBMethod(double eps, const std::function<double(double)> &f) : IntegrationMethod(eps, f) {}

double MPI_TBBMethod::tbb_run(double a, double b) {
	double x = simpson(a, b, 4);
    double y = simpson(a, b, 8);
    if (fabs(x - y) < eps) {
        return y;
    } else {
        tbb::task_group group;

        group.run([&] {
            x = tbb_run(a, (a + b) / 2.0);
        });

        group.run_and_wait([&] {
            y = tbb_run((a + b) / 2.0, b);
        });

        return x + y;
    }
}

double MPI_TBBMethod::integrate(double a, double b) {
	double answer = 0.0;
	int nodes = 1;
    MPI_Comm_size(MPI_COMM_WORLD, &nodes);

    int node_id = 0;
    MPI_Comm_rank(MPI_COMM_WORLD, &node_id);

	double h = (b - a) / (double) nodes;

	a = a + h * node_id;
    b = a + h;

    double node_answer = tbb_run(a, b);
    MPI_Barrier(MPI_COMM_WORLD);
    if (node_id == 0) {
    	answer += node_answer;
        for (int i = 1; i < nodes; i++) {
            MPI_Recv(&node_answer, 1, MPI_DOUBLE, MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
            answer += node_answer;
        }
    } else {
        MPI_Send(&node_answer, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
    }

    MPI_Barrier(MPI_COMM_WORLD);

	return answer;
}

MPI_TBBMethod::~MPI_TBBMethod() {}

В данной реализации отрезок [a, b] делится на число доступных процессоров, каждый отрезок считается на своем процессоре, рекурсия заменяется на систему задач библиотеки Intel MPI для более эффективного использования многоядерных процессоров.

3.3.2 Реализация только с помощью MPI

Исходный код

В данной реализации рекурсия заменяется на систему задач, каждая задача представляет собой подсчет интеграла на отрезке, в случае если необходимая точность не достигнута (см. Математическое описание алгоритма), создаются две новые задачи с меньшими отрезками и распределяются между другими процессорами. Данная реализация в теории решает проблему неравномерности гладких и не гладких участков функции, в случае когда большая часть участков гладкая и лишь небольшой участок функции сильно осциллирует.

3.4 Тесты производительности

Все тесты проводились на суперкомпьютере Ломоносов.
Тестировались две функции разной степени сложности.
Из тестов видно, что параллельная реализация с использованием технологий MPI и Intel TBB показывает прекрасную масштабируемость, в то же время версия на чистом MPI с управляющим узлом сильно проигрывает. Это связано с тем, что при синхронизации задач между процессорами происходит много пересылок сообщений MPI, в сравнении с вычислением значений функции в узловых точках. Это сильно замедляет работу программы в целом, более того при увеличении количества узлов управляющий узел не справляется с очень большим потоком сообщений от вычислительных узлов.

3.4.1 MPI + TBB

Рисунок 1. Результаты масштабируемости для первой функции

Рисунок 2. Результаты масштабируемости для второй функции

3.4.2 MPI Farm

Рисунок 3. Результаты масштабируемости для первой функции

Рисунок 4. Результаты масштабируемости для второй функции