Участник:Joinmek/Вычисление определенного интеграла с использованием адаптивно сгущающейся сетки

Вычисление определенного интеграла с использованием адаптивно сгущающейся сетки
Последовательный алгоритм
Объём входных данных	[math]O(1)[/math]
Объём выходных данных	[math]O(1)[/math]

Основные авторы описания: Можарова В. А.

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

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

Основная идея большинства методов численного интегрирования состоит в замене подынтегральной функции на более простую, интеграл от которой легко вычисляется аналитически. При замене подынтегральной функции на полином нулевой, первой и второй степени получаются соответственно методы прямоугольников, трапеций и парабол (Симпсона).

Рис.1. Метод трапеций

Чтобы посчитать определенный интеграл, участок интегрирования делится на отрезки, на которых подсчитывается значение интеграла от новой (упрощенной) функции, после чего значения интегралов суммируются на всех отрезках. Но чтобы посчитать интеграл с определенной точностью, мы не можем заранее узнать, какие размеры отрезков, нужно использовать, т. к. это зависит от самой функции. Для этого используют адаптивное разбиение участка интегрирования: размер конкретного отрезка определяется во время подсчета значения функции на нем самом. Если заданная точность не была достигнута на данном отрезке, он рекурсивно разбивается на несколько новых и функция подсчитывается уже на них. Таким образом, выигрыш от применения рассмотренного алгоритма рекурсивного разбиения достигается за счет возможности использования сеток с разным числом узлов на разных участках интервала интегрирования.

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

Пусть у нас есть интеграл, который надо вычислить с точностью ε:

Предположим, что в результате адаптации сетки для данной функции, на отрезке [A, B] получилась сетка, содержащая n + 1 узел. Тогда согласно методу трапеций, можно численно найти определенный интеграл от функции f(x) на отрезке [A, B]:

Будем полагать значение J найденным с точностью ε, если выполнено условие:

где n₁ и n₂ - количество узлов, образующих две разные сетки. Пользуясь формулой

рекурсивно делим отрезки интегрирования до достижения на них заданной точности.

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

Основное время работы алгоритма приходится на вычисление площадей на отрезках, получившихся в результате рекурсивного разбиения:

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

Алгоритм итеративный. Сначала значение текущего отрезка равно общему отрезу интегрирования. На каждой итерации:

1. На текущем отрезке считается приближенное значение интеграла

2. Текущий отрезок делится пополам, и на получившихся отрезках высчитываются новые интегралы

3. Если достигнута заданная точность, то

3.1 Добавить значение интеграла текущего отрезка к общему интегралу

3.2 Если стек пуст, то завершить алгоритм

3.3 Достать из стека отрезок и присвоить его значение текущему отрезку.

4. Иначе

4.1 Один из получившихся маленьких отрезков положить в стек

4.2 Текущему отрезку присвоить значение второго маленького отрезка

1.5 Схема реализации алгоритма

Для хранения отрезков интегрирования в данном алгоритме понадобится стек .

   sp=0  //    указатель вершины стека
   struct {
         A, B, fA, fB, s
   } stk[1000] 

У стека определен следующий набор операций:

1. PUT_INTO_STACK - положить элемент в стек

   #define PUT_INTO_STACK(A, B, fA, fB, s)
   {
       stk[sp].A=A
       stk[sp].B=B
       stk[sp].fA=fA
       stk[sp].fB=fB
       stk[sp].s=s
       sp++
   } 

2. GET_FROM_STACK - взять верхний элемент из стека

   #define GET_FROM_STACK(A, B, fA, fB, s)
   {
       sp--
       A=stk[sp].A
       B=stk[sp].B
       fA=stk[sp].fA
       fB=stk[sp].fB
       s=stk[sp].s
   }

3. STACK_IS_NOT_FREE - проверить стек на наличие в нем элементов

   #define STACK_IS_NOT_FREE (sp > 0)

   LocalStack (A, B)
   {
       // В переменной J храним ответ
       J=0
       // Вычисляем значение функции на границах интегрируемого отрезка
       fA = f(A)
       fB = f(B)
       // Вычисляем площадь с помощью метода трапеций
       sAB = (fA + fB) * (B - A) / 2
       while (1)
       {
           C = (A + B) / 2  // Находим середину отрезка
           fc = f(C)  // Вычисляем в середине значение функции
           // Считаем приближенное значение интеграла на получившихся отрезках
           sAC = (fA + fC) * (C - A) / 2
           sCB = (fC + fB ) * (B-C) / 2
           sACB = sAC + sCB
           if (|sAB - sACB| ≥ ε|sAB|) {  // Если заданная точность не достигнута
               PUT_INTO_STACK(A, C, fA, fC, sAC)  // то один из отрезков кладем в стек
               A = C  // и решаем задачу для оставшегося отрезка
               fA = fC
               sAB = sCB
           } else  // Если заданная точность достигнута
           {
               J += sACB  // то к общей сумме добавляем полученное число
               if (STACK_IS_NOT_FREE)  // если стек пуст, то цикл завершается - ответ получен
                   break 
               GET_FROM_STACK(A, B, fA, fB, sAB)  // иначе берем из стека новый отрезок
           } 
       }
       return J 
   }

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

Количество делений отрезков интегрирования зависит от сложности входной функции, поэтому нельзя оценить сложность данного алгоритма. Если предположить, что алгоритм будет останавливаться при достижении длины отрезка меньше некоторого порога G, то сложность можно оценить как O((B - A) / G).

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

В параллельной версии алгоритма, появляется дополнительный глобальный стек, доступ к которому имеют все процессы. Он находится в общей памяти и так же содержит границы отрезков. Процессы обмениваются с ним данными, распределяя таким образом отрезки между собой. Когда в локальном стеке заканчиваются элементы, процесс обращается к глобальному стеку, чтобы получить их. Причем любой процесс может отдавать ему отрезки, если глобальный стек пуст в то время, как локальный имеет более одного отрезка. Изначально в глобальный стек помещается один отрезок, соответствующий участку интегрирования. Существуют реализации, где локальные стеки полностью отсутствуют, а все отрезки хранятся в глобальном стеке. Такой вариант проще реализовать, но его недостаток в том, что увеличивается количество обращений процессов к глобальному стеку, и каждый процесс должен ждать своей очереди, что может привести к увеличению времени работы программы.

Рис.2. Информационный граф для алгоритма глобального стека. Ch - выбор отрезка интегрирования, Sq - вычисление площади на отрезке, Div - деление отрезка пополам, Tol - проверка достижения точности, Wr - запись отрезка в лок./глоб. стек, Sum - суммирование результатов

1.8 Ресурс параллелизма алгоритма

Сложность параллельного алгоритма нельзя оценить из-за зависимости от сложности интегрируемой функции так же, как и у последовательного алгоритма. Но если наложить ограничение на минимально допустимую величину отрезка, как было сделано ранее, то его можно оценить следующим образом: O((B - A) / (G * p)), где p -количество процессов. Каждый процесс считает интеграл на своем отрезке вне зависимости от других процессов, влияние других процессоров происходит исключительно в момент ожидания получения доступа к глобальному стеку.

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

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

1.10 Свойства алгоритма

1. Алгоритм адаптируется под сложность функции

2. Для каждого узла сетки значение функции высчитывается один раз

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

4. Алгоритм является детерминированным. На это свойство влияет выполнение свойства ассоциативности сложения, но не влияют свойства интегрируемой функции.

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

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

7. Данный алгоритм является устойчивым не для всех функций. Например, для sinx алгоритм неустойчив. Можно рассмотреть отрезки интегрирования [0, 2*pi] и [0, 2*pi + delta]. В первом случае, алгоритм остановится на первой итерации, и значение интеграла будет равно 0. При изменении правой границы отрезка на малое delta, алгоритм не остановится на первой итерации, и значение интеграла будет сильно отличаться от 0.

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

2.4 Масштабируемость алгоритма и его реализации

Автором данной статьи была выбрана реализация RWTH AACHEN UNIVERSITY на языке C.

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

Чтобы протестировать систему была выбрана функция sin(x). Вычисление интеграла запускается с разными значениями точности, чтобы пронаблюдать за изменением времени в зависимости от количества итоговых отрезков. Программа запускалась на Blue Gene. Использовался компилятор mpixlcxx_r с опцией -O2. Внутри программы использовалась стандартная библиотека языка C - math, которая была выбрана для подсчета функции синуса в заданных точках, а также библиотека OpenMPI версии 1.68.

Рис.3. Зависимость времени работы алгоритма от числа процессов и заданной точности

Рис.4. Зависимость эффективности алгоритма от числа процессов и заданной точности

Масштабируемость:

Система запускалась на количестве процессов от 8 до 128, начиная с 16 процессов шаг был равен 16. Точность градировалась от 0.00001 до 0.0001 с величиной шага 0.00001. Максимальная эффективность была достигнута на 8 процессах и равна 0.9389. Минимальная - на 128, ее величина 0.2914.

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

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

Для оценки масштабируемости алгоритма на системе с общей памятью, реализация RWTH AACHEN UNIVERSITY была запущена на архитектуре IBM pSeries 690 HPC Regatta. Программа запускалась на 1, 2, 4, 8 и 16 потоках. Точность так же, как и в первом случае, градировалась от 0.00001 до 0.0001 с величиной шага 0.00001. Максимальная эффективность была достигнута 2 потоках и равна 0.9504; минимальная - на 16, ее величина - 0.67484. При увеличении количества потоков эффективность падает из-за конкуренции доступа к глобальному стеку. При увеличении сложности для одинакового числа процессов эффективность остается примерно на том же уровне.

Рис.3. Зависимость времени работы алгоритма от числа потоков и заданной точности

Рис.3. Зависимость времени работы алгоритма от числа процессов и заданной точности

2.7 Существующие реализации алгоритма

В России данный алгоритм реализован Якобовским М. В.

За границей осуществила реализацию компания RWTH Compute Cluster; группа исследователей во главе с Kamesh Arumugam; а также команда Berntsen, J.; Espelid, T. O.; and Genz, A..

3 Литература

<references \>