Адаптивное интегрирование

Незаконченный вариант статьи

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

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

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

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

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

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

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

1.3 Макростуктура алгоритма

Основную часть метода составляют рекурсивные вызовы интегралов отрезков меньшей длинны.

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

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

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

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

2 Программная реализация алгоритма

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

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

2.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() {}

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

2.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 для более эффективного использования многоядерных процессоров.

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

Исходный код

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

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

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

2.4.1 MPI + TBB

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

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

2.4.2 MPI Farm

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

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